TABLE OF CONTENTS


ABINIT/m_ddb [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddb

FUNCTION

  This module contains the declaration of data types and methods
  used to handle the blocks of data in DDB files:
  blkval, nrm, qpt, flg, and associated dimensions
  Main entry point for client code that needs to read the DDB data.

COPYRIGHT

 Copyright (C) 2011-2018 ABINIT group (MJV, XG, MT, MM, MVeithen, MG, PB, JCC, SP)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_ddb
29 
30  use defs_basis
31  use defs_abitypes
32  use defs_datatypes
33  use m_abicore
34  use m_errors
35  use m_xmpi
36  use m_ddb_hdr
37  use m_dtset
38 
39  use m_fstrings,       only : sjoin, itoa, ktoa
40  use m_numeric_tools,  only : mkherm
41  use m_symtk,          only : mati3inv, matr3inv, littlegroup_q, symatm
42  use m_io_tools,       only : get_unit
43  use m_copy,           only : alloc_copy
44  use m_geometry,       only : phdispl_cart2red, mkrdim, xred2xcart, metric
45  use m_crystal,        only : crystal_t, crystal_init
46  use m_dynmat,         only : cart29, d2sym3, cart39, d3sym, chneu9, asria_calc, asria_corr, asrprs, &
47 &                             dfpt_phfrq, sytens
48  use m_pawtab,         only : pawtab_type, pawtab_nullify, pawtab_free
49  use m_psps,           only : psps_copy, psps_free
50 
51  implicit none
52 
53  private
54 
55  public :: ddb_write_blok   ! Writes blocks of data in the DDBs.
56  public :: dfptnl_doutput   ! Write the matrix of third-order derivatives to the output file and the DDB
57  public :: gtblk9           ! Finds the block containing the derivatives of the total energy.
58  public :: read_blok8       ! This routine reads blocks of data in the DDBs.
59  public :: rdddb9           ! This routine reads the derivative database entirely,
60  public :: nlopt            ! Output of all quantities related to third-order derivatives of the energy.
61  public :: chkin9
62  public :: carttransf       ! Transform a second-derivative matrix (EIG2D) from reduced
63                             ! coordinates to cartesian coordinates.
64 
65  integer,public,parameter :: DDB_VERSION=100401
66  ! DDB Version number.
67  ! TODO: Remove other occurrences of this magic number.
68  ! Postponed to avoid conflicts with Samuel's branch in dfpt_looppert
69 
70  real(dp),public,parameter :: DDB_QTOL=2.0d-8
71  ! Tolerance for the identification of two wavevectors

m_db_blk/read_blok8 [ Functions ]

[ Top ] [ Functions ]

NAME

 read_blok8

FUNCTION

 This routine reads blocks of data in the DDBs.

INPUTS

 mpert =maximum number of ipert
 msize=maximum size of the arrays flags and values
 nunit=unit number for the data block file

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 ddb = ddb block datastructure
 ddb%typ=type of the block:
   0 => total energy
   1 => second-order energy derivatives, non-stationary block
   2 => second-order energy derivatives, stationary block
   3 => third-order energy derivatives
   4 => first-order energy derivatives: forces, stresses and polarization
   5 => second-order eigenvalue derivatives
 ddb%flg(msize)=flag for every matrix element (0=> the element is
  not in the data block), (1=> the element is in the data blok)
 ddb%qpt(9)=wavevector of the perturbation(s). The elements from
  1 to 3 are used if we are dealing with the 2nd derivative of
  total energy (only one wavevector), while all elements are
  used in case of a third order derivative of total energy (three wavevector could be present)
 ddb%nrm(3)=normalization factors for the three allowed wavevectors.
 ddb%val(2,msize)=real(dp), complex, value of the matrix elements that are present in the data block
 [blkval2(2,msize,mband,nkpt)]= value of the matrix elements that are present in a block of EIGR2D/EIGI2D

NOTES

 only executed by one processor.

PARENTS

      m_ddb,mblktyp1,mblktyp5,thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

 877 subroutine read_blok8(ddb,iblok,mband,mpert,msize,nkpt,nunit,&
 878 &     blkval2,kpt) !optional
 879 
 880 
 881 !This section has been created automatically by the script Abilint (TD).
 882 !Do not modify the following lines by hand.
 883 #undef ABI_FUNC
 884 #define ABI_FUNC 'read_blok8'
 885 !End of the abilint section
 886 
 887  implicit none
 888 
 889 !Arguments -------------------------------
 890 !scalars
 891  integer,intent(in) :: mband,mpert,msize,nkpt,nunit
 892  integer, intent(in) :: iblok
 893  type(ddb_type),intent(inout) :: ddb
 894 !arrays
 895  real(dp),intent(out),optional :: kpt(3,nkpt)
 896  real(dp),intent(out),optional :: blkval2(2,msize,mband,nkpt)
 897 
 898 !Local variables -------------------------
 899 !scalars
 900  integer :: band,iband,idir1,idir2,idir3,ii,ikpt,index,ipert1,ipert2,ipert3,nelmts
 901  real(dp) :: ai,ar
 902  character(len=32) :: name
 903  character(len=500) :: message
 904 
 905 ! *********************************************************************
 906 
 907 !Zero every flag
 908  ddb%flg(1:msize, iblok)=0
 909  if(present(blkval2))blkval2(:,:,:,:)=zero
 910  if(present(kpt))kpt(:,:)=zero
 911 
 912 !Read the block type and number of elements
 913  read(nunit,*)
 914  read(nunit, '(a32,12x,i8)' )name,nelmts
 915 
 916  if(name==' 2nd derivatives (non-stat.)  - ' .or. name==' 2rd derivatives (non-stat.)  - ')then
 917    ddb%typ(iblok)=1
 918  else if(name==' 2nd derivatives (stationary) - ' .or. name==' 2rd derivatives (stationary) - ')then
 919    ddb%typ(iblok)=2
 920  else if(name==' 3rd derivatives              - ')then
 921    ddb%typ(iblok)=3
 922  else if(name==' Total energy                 - ')then
 923    ddb%typ(iblok)=0
 924  else if(name==' 1st derivatives              - ')then
 925    ddb%typ(iblok)=4
 926  else if(name==' 2nd eigenvalue derivatives   - ' .or. name==' 2rd eigenvalue derivatives   - ')then
 927    ddb%typ(iblok)=5
 928  else
 929    write(message,'(6a)')&
 930    'The following string appears in the DDB in place of',&
 931    ' the block type description :',ch10,trim(name),ch10,&
 932    'Action: check your DDB.'
 933    MSG_ERROR(message)
 934  end if
 935 
 936 !Read the 2nd derivative block
 937  if(ddb%typ(iblok)==1.or.ddb%typ(iblok)==2)then
 938 
 939 !  First check if there is enough space to read it
 940    if(msize<(3*mpert*3*mpert))then
 941      write(message,'(3a)')&
 942      'There is not enough space to read a second-derivative block.',ch10,&
 943      'Action: increase msize and recompile.'
 944      MSG_ERROR(message)
 945    end if
 946 
 947 !  Read the phonon wavevector
 948    read(nunit, '(4x,3es16.8,f6.1)' )(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
 949 
 950 !  Read every element
 951    do ii=1,nelmts
 952      read(nunit,*)idir1,ipert1,idir2,ipert2,ar,ai
 953      index=idir1+3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
 954      ddb%flg(index,iblok)=1
 955      ddb%val(1,index,iblok)=ar
 956      ddb%val(2,index,iblok)=ai
 957    end do
 958 
 959 !  Read the 3rd derivative block
 960  else if(ddb%typ(iblok)==3)then
 961 
 962 !  First check if there is enough space to read it
 963    if(msize<(3*mpert*3*mpert*3*mpert))then
 964      write(message, '(a,a,a,i10,a,i10,a,a,a)' )&
 965      'There is not enough space to read a third-derivative block.',ch10,&
 966      'The size provided is only ',msize,' although ',3*mpert*3*mpert*3*mpert,' is needed.',ch10,&
 967      'Action: increase msize and recompile.'
 968      MSG_ERROR(message)
 969    end if
 970 
 971 !  Read the perturbation wavevectors
 972    read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
 973    read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=4,6),ddb%nrm(2,iblok)
 974    read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=7,9),ddb%nrm(3,iblok)
 975 
 976 !  Read every element
 977    do ii=1,nelmts
 978      read(nunit,'(6i4,2d22.14)')idir1,ipert1,idir2,ipert2,idir3,ipert3,ar,ai
 979      index=idir1+                     &
 980 &     3*((ipert1-1)+mpert*((idir2-1)+ &
 981 &     3*((ipert2-1)+mpert*((idir3-1)+3*(ipert3-1)))))
 982      ddb%flg(index,iblok)=1
 983      ddb%val(1,index,iblok)=ar
 984      ddb%val(2,index,iblok)=ai
 985    end do
 986 
 987 !  Read the total energy
 988  else if(ddb%typ(iblok)==0)then
 989 
 990 !  First check if there is enough space to read it
 991    if(msize<1)then
 992      write(message, '(3a,i0,3a)' )&
 993 &     'There is not enough space to read a total energy block.',ch10,&
 994 &     'The size provided is only ',msize,' although 1 is needed.',ch10,&
 995 &     'Action: increase msize and recompile.'
 996      MSG_ERROR(message)
 997    end if
 998 
 999 !  Read the total energy
1000    read(nunit,'(2d22.14)')ar,ai
1001    ddb%flg(1,iblok)=1
1002    ddb%val(1,1,iblok)=ar
1003    ddb%val(2,1,iblok)=ai
1004 
1005 !  Read the 1st derivative block
1006  else if (ddb%typ(iblok) == 4) then
1007 
1008 !  First check if there is enough space to read it
1009    if (msize < (3*mpert)) then
1010      write(message, '(3a,i0,a,i0,3a)' )&
1011      'There is not enough space to read a first-derivative block.',ch10,&
1012      'The size provided is only ',msize,' although ',3*mpert,' is needed.',ch10,&
1013      'Action: increase msize and recompile.'
1014      MSG_ERROR(message)
1015    end if
1016 
1017 !  Read every element
1018    do ii=1,nelmts
1019      read(nunit,'(2i4,2d22.14)')idir1,ipert1,ar,ai
1020      index=idir1 + 3*(ipert1 - 1)
1021      ddb%flg(index,iblok)=1
1022      ddb%val(1,index,iblok)=ar
1023      ddb%val(2,index,iblok)=ai
1024    end do
1025 
1026 !  Read the 2nd eigenvalue derivative block
1027  else if(ddb%typ(iblok)==5)then
1028 
1029 !  First check if there is enough space to read it
1030    if(msize<(3*mpert*3*mpert))then
1031      write(message, '(3a,i0,a,i0,3a)' )&
1032      'There is not enough space to read a second-derivative block.',ch10,&
1033      'The size provided is only ',msize,' although ',3*mpert*3*mpert*mband*nkpt,' is needed.',ch10,&
1034      'Action: increase msize and recompile.'
1035      MSG_ERROR(message)
1036    end if
1037 
1038 !  Read the phonon wavevector
1039    read(nunit, '(4x,3es16.8,f6.1)' )(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
1040 
1041 !  Read the K point and band
1042    if(present(blkval2).and.present(kpt))then
1043      do ikpt=1,nkpt
1044        read(nunit, '(9x,3es16.8)')(kpt(ii,ikpt),ii=1,3)
1045        do iband=1,mband
1046          read(nunit, '(6x,i3)') band
1047 !        Read every element
1048          do ii=1,nelmts
1049            read(nunit,*)idir1,ipert1,idir2,ipert2,ar,ai
1050            index=idir1+3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1051            ddb%flg(index,iblok)=1
1052            blkval2(1,index,iband,ikpt)=ar
1053            blkval2(2,index,iband,ikpt)=ai
1054          end do !nelmts
1055        end do  !band
1056      end do !kpt
1057    end if
1058  end if
1059 
1060 end subroutine read_blok8

m_ddb/asr_t [ Types ]

[ Top ] [ m_ddb ] [ Types ]

NAME

  asr_t

FUNCTION

  Object used to enforce the acoustic sum rule from the Dynamical matrix at Gamma.
  Wraps several approaches that can be activated via the `asr` option.

SOURCE

173  type,public :: asrq0_t
174 
175    integer :: iblok = 0
176     ! Index of the Gamma block in the DDB.
177     ! Set to 0 if no block was found. Client code can use this flag to understand
178     ! if ASR can be enforced.
179 
180    integer :: asr
181    ! Option for the application of the ASR (input variable).
182 
183    integer :: natom
184     ! Number of atoms.
185 
186    real(dp),allocatable :: d2asr(:,:,:,:,:)
187    ! d2asr,(2,3,natom,3,natom))
188    ! In case the interatomic forces are not calculated, the
189    ! ASR-correction (d2asr) has to be determined here from the Dynamical matrix at Gamma.
190 
191    ! singular, uinvers and vtinvers are allocated and used only if asr in [3,4]
192    ! i.e. Rotational invariance for 1D and 0D systems. dims=3*natom*(3*natom-1)/2
193    real(dp),allocatable :: singular(:)
194    ! singular,(1:dims))
195 
196    real(dp),allocatable :: uinvers(:,:)
197    ! uinvers,(1:dims,1:dims))
198 
199    real(dp),allocatable :: vtinvers(:,:)
200    ! vtinvers,(1:dims,1:dims))
201 
202  end type asrq0_t
203 
204  public :: asrq0_apply      ! Impose the acoustic sum rule based on the q=0 block found in the DDB file.
205  public :: asrq0_free       ! Free memory

m_ddb/asrq0_apply [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 asrq0_apply

FUNCTION

  Impose the acoustic sum rule based on the q=0 block found in the DDB file.

INPUTS

  asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
  natom=Number of atoms per unit cell.
  mpert=Maximum number of perturbation (reported in ddb%mpert)
  msize=Maximum size of array ddb%val
  xcart(3,natom)=Atomic positions in Cartesian coordinates

SIDE EFFECTS

   d2cart=matrix of second derivatives of total energy, in cartesian coordinates
   Input: Values stored in ddb%
   Output: Changed to enforce ASR.

PARENTS

      anaddb,m_ddb,m_phonons

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

2957 subroutine asrq0_apply(asrq0, natom, mpert, msize, xcart, d2cart)
2958 
2959 
2960 !This section has been created automatically by the script Abilint (TD).
2961 !Do not modify the following lines by hand.
2962 #undef ABI_FUNC
2963 #define ABI_FUNC 'asrq0_apply'
2964 !End of the abilint section
2965 
2966  implicit none
2967 
2968 !Arguments ------------------------------------
2969 !scalars
2970  integer,intent(in) :: natom, msize, mpert
2971  type(asrq0_t),intent(inout) :: asrq0
2972 !arrays
2973  real(dp),intent(in) :: xcart(3,natom)
2974  real(dp),intent(inout) :: d2cart(2,msize)
2975 
2976 ! ************************************************************************
2977 
2978  if (asrq0%asr /= 0 .and. asrq0%iblok == 0) then
2979    MSG_WARNING("asr != 0 but DDB file does not contain q=Gamma. D(q) cannot be corrected")
2980    return
2981  end if
2982 
2983  select case (asrq0%asr)
2984  case (0)
2985    return
2986  case (1,2,5)
2987    call asria_corr(asrq0%asr, asrq0%d2asr, d2cart, mpert, natom)
2988  case (3,4)
2989    ! Impose acoustic sum rule plus rotational symmetry for 0D and 1D systems
2990    call asrprs(asrq0%asr,2,3,asrq0%uinvers,asrq0%vtinvers,asrq0%singular,d2cart,mpert,natom,xcart)
2991  case default
2992    MSG_ERROR(sjoin("Wrong value for asr:", itoa(asrq0%asr)))
2993  end select
2994 
2995 end subroutine asrq0_apply

m_ddb/asrq0_free [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 asrq0_free

FUNCTION

   Free dynamic memory

PARENTS

      anaddb,m_effective_potential_file

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

3015 subroutine asrq0_free(asrq0)
3016 
3017 
3018 !This section has been created automatically by the script Abilint (TD).
3019 !Do not modify the following lines by hand.
3020 #undef ABI_FUNC
3021 #define ABI_FUNC 'asrq0_free'
3022 !End of the abilint section
3023 
3024  implicit none
3025 
3026 !Arguments ------------------------------------
3027  type(asrq0_t),intent(inout) :: asrq0
3028 
3029 ! ************************************************************************
3030 
3031  ! real
3032  if (allocated(asrq0%d2asr)) then
3033    ABI_FREE(asrq0%d2asr)
3034  end if
3035 
3036  if (allocated(asrq0%singular)) then
3037    ABI_FREE(asrq0%singular)
3038  end if
3039 
3040  if (allocated(asrq0%uinvers)) then
3041    ABI_FREE(asrq0%uinvers)
3042  end if
3043 
3044  if (allocated(asrq0%vtinvers)) then
3045    ABI_FREE(asrq0%vtinvers)
3046  end if
3047 
3048 end subroutine asrq0_free

m_ddb/carteig2d [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 carteig2d

FUNCTION

 Transform a second-derivative matrix (EIG2D) from reduced
 coordinates to cartesian coordinates

INPUTS

  blkflg(3,mpert,3,mpert,nblok)=
   ( 1 if the element of the dynamical matrix has been calculated ;
     0 otherwise )
  blkval(2,3,mpert,3,mpert,nblok)=DDB values
  gprimd(3,3)=basis vector in the reciprocal space
  iblok=number of the blok that will be transformed
  mpert =maximum number of ipert
  natom=number of atom
  nblok=number of blocks (dimension of blkflg and blkval)
  rprimd(3,3)=basis vector in the real space

OUTPUT

  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
  2DTE matrix has been calculated correctly ; 0 otherwise )
  d2cart(2,3,mpert,3,mpert)=
    dynamical matrix, effective charges, dielectric tensor,....
    all in cartesian coordinates

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1979 subroutine carteig2d(blkflg,blkval,carflg,d2cart,&
1980 & gprimd,iblok,mpert,natom,nblok,rprimd)
1981 
1982 
1983 !This section has been created automatically by the script Abilint (TD).
1984 !Do not modify the following lines by hand.
1985 #undef ABI_FUNC
1986 #define ABI_FUNC 'carteig2d'
1987 !End of the abilint section
1988 
1989  implicit none
1990 
1991 !Arguments -------------------------------
1992 !scalars
1993  integer,intent(in) :: iblok,mpert,natom,nblok
1994 !arrays
1995  integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok)
1996  integer,intent(out) :: carflg(3,mpert,3,mpert)
1997  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3)
1998  real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert)
1999 
2000 !Local variables -------------------------
2001 !scalars
2002  integer :: idir1,idir2,ii,ipert1,ipert2
2003 !arrays
2004  integer :: flg1(3),flg2(3)
2005  real(dp) :: vec1(3),vec2(3)
2006 
2007 ! *********************************************************************
2008 
2009 !First, copy the data blok in place.
2010  d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok)
2011 
2012 !Cartesian coordinates transformation (in two steps)
2013 !First step
2014  do ipert1=1,mpert
2015    do ipert2=1,mpert
2016      do ii=1,2
2017        do idir1=1,3
2018          do idir2=1,3
2019            vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2)
2020 !          Note here blkflg
2021            flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok)
2022          end do
2023          call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2)
2024          do idir2=1,3
2025            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2)
2026 !          And here carflg
2027            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2)
2028          end do
2029        end do
2030      end do
2031    end do
2032  end do
2033 
2034 !Second step
2035  do ipert1=1,mpert
2036    do ipert2=1,mpert
2037      do ii=1,2
2038        do idir2=1,3
2039          do idir1=1,3
2040            vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2)
2041 !          Note here carflg
2042            flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2)
2043          end do
2044          call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2)
2045          do idir1=1,3
2046            d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1)
2047 !          And here carflg again
2048            carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1)
2049          end do
2050        end do
2051      end do
2052    end do
2053  end do
2054 
2055 end subroutine carteig2d

m_ddb/carttransf [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 carttransf

FUNCTION

 Transform a second-derivative matrix (EIG2D) from reduced
 coordinates to cartesian coordinates.

INPUTS

  blkflg(msize,nblok)=
   ( 1 if the element of the dynamical matrix has been calculated ;
     0 otherwise )
  gprimd(3,3)=basis vector in the reciprocal space
  iqpt  = number of the Q-point currently used
  mband = maximal number of bands
  mpert = maximum number of ipert
  msize = size of the EIG2D arrays (3*mpert*3*mpert)
  natom = number of atom
  nblok = number of bloks in blkflg
  nkpt  = number of K-points
  rprimd(3,3) = basis vector in the real space

OUTPUT

  carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
  EIG2D matrix has been calculated correctly ; 0 otherwise )

 SIDE EFFECT
 blkval2(2,msize,mband,nkpt)=Second order eigenvalues (EIG2D)
 is transformed from reduced coordinates to cartesian coordinates

PARENTS

      thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1859 subroutine carttransf(blkflg,blkval2,carflg,gprimd,iqpt,mband,&
1860 & mpert,msize,natom,nblok,nkpt,rprimd)
1861 
1862 
1863 !This section has been created automatically by the script Abilint (TD).
1864 !Do not modify the following lines by hand.
1865 #undef ABI_FUNC
1866 #define ABI_FUNC 'carttransf'
1867 !End of the abilint section
1868 
1869  implicit none
1870 
1871 !Arguments ------------------------------------
1872 !scalars
1873  integer,intent(in) :: mband,msize
1874  integer,intent(in) :: iqpt
1875  integer,intent(in) :: mpert,nblok
1876  integer,intent(inout) :: natom,nkpt
1877 !arrays
1878  integer,intent(in) :: blkflg(msize,nblok)
1879  integer,intent(out) :: carflg(3,mpert,3,mpert)
1880  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
1881  real(dp),intent(inout) :: blkval2(2,msize,mband,nkpt)
1882 
1883 !Local variables-------------------------------
1884 !scalars
1885 integer :: iatom1,iatom2,iband,idir1,idir2,ikpt
1886 integer :: index
1887 !arrays
1888 real(dp),allocatable :: blkflgtmp(:,:,:,:,:)
1889 real(dp),allocatable :: blkval2tmp(:,:,:,:,:,:)
1890 real(dp),allocatable :: d2cart(:,:,:,:,:)
1891 
1892 ! *********************************************************************
1893 
1894 !Start by allocating local arrays
1895  ABI_MALLOC(blkflgtmp,(3,mpert,3,mpert,1))
1896  ABI_MALLOC(blkval2tmp,(2,3,mpert,3,mpert,1))
1897  ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
1898 
1899 !Begin by formating the arrays to be compatible with cart29
1900 !Then call cart29 to transform the arrays in cartesian coordinates
1901 !Finally reformat the cartesian arrays in old format
1902  do ikpt=1,nkpt
1903    do iband=1,mband
1904 
1905      do idir1=1,3
1906        do iatom1=1,mpert
1907          do idir2=1,3
1908            do iatom2=1,mpert
1909              index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1)))
1910              blkflgtmp(idir1,iatom1,idir2,iatom2,1) = blkflg(index,iqpt)
1911              blkval2tmp(:,idir1,iatom1,idir2,iatom2,1) = blkval2(:,index,iband,ikpt)
1912            end do
1913          end do
1914        end do
1915      end do
1916 
1917 !    The 1sin the argument of cart29 are respectively iblok and nblok. We are doing only one blok.
1918      call carteig2d(blkflg(:,iqpt),blkval2tmp,carflg,d2cart,gprimd,1,mpert,natom,1,rprimd)
1919 
1920      do idir1=1,3
1921        do iatom1=1,mpert
1922          do idir2=1,3
1923            do iatom2=1,mpert
1924              index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1)))
1925              blkval2(:,index,iband,ikpt) = d2cart(:,idir1,iatom1,idir2,iatom2)
1926            end do
1927          end do
1928        end do
1929      end do
1930 
1931    end do
1932  end do
1933 
1934 !Deallocating local arrays
1935  ABI_FREE(blkflgtmp)
1936  ABI_FREE(blkval2tmp)
1937  ABI_FREE(d2cart)
1938 
1939 end subroutine carttransf

m_ddb/chkin9 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 chkin9

FUNCTION

 Check the value of some input parameters.
 Send error message and stop if needed.
 Also transform the meaning of atifc

INPUTS

 atifc(natom)=list of the atom ifc to be analysed
 natifc= number of atom ifc to be analysed
 natom= number of atoms

OUTPUT

 atifc(natom) =  atifc(ia) equals 1 if the analysis of ifc
  has to be done for atom ia; otherwise 0.

NOTES

 Only for one processor (no use of wrtout)

PARENTS

      m_ddb,thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1380 subroutine chkin9(atifc,natifc,natom)
1381 
1382 
1383 !This section has been created automatically by the script Abilint (TD).
1384 !Do not modify the following lines by hand.
1385 #undef ABI_FUNC
1386 #define ABI_FUNC 'chkin9'
1387 !End of the abilint section
1388 
1389  implicit none
1390 
1391 !Arguments -------------------------------
1392 !scalars
1393  integer,intent(in) :: natifc,natom
1394 !arrays
1395  integer,intent(inout) :: atifc(natom)
1396 
1397 !Local variables -------------------------
1398 !scalars
1399  integer :: iatifc
1400  character(len=500) :: message
1401 !arrays
1402  integer,allocatable :: work(:)
1403 
1404 ! *********************************************************************
1405 
1406  if(natifc>natom)then
1407    write(message, '(a,i0,3a,i0,3a)' )&
1408 &   'The number of atom ifc in the input files',natifc,',',ch10,&
1409 &   'is larger than the number of atoms',natom,'.',ch10,&
1410 &   'Action: change natifc in the input file.'
1411    MSG_ERROR(message)
1412  end if
1413 
1414  if(natifc>=1)then
1415    ABI_MALLOC(work,(natom))
1416    work(:)=0
1417 
1418    do iatifc=1,natifc
1419      if(atifc(iatifc)<=0.or.atifc(iatifc)>natom)then
1420        write(message, '(a,i0,5a,i0,3a)' )&
1421 &       'For iatifc=',iatifc,', the number of the atom ifc to be ',ch10,&
1422 &       'analysed is not valid : either negative, ',ch10,&
1423 &       'zero, or larger than natom =',natom,'.',ch10,&
1424 &       'Action: change atifc in your input file.'
1425        MSG_ERROR(message)
1426      end if
1427      work(atifc(iatifc))=1
1428    end do
1429 
1430    atifc(1:natom)=work(:)
1431    ABI_FREE(work)
1432  end if
1433 
1434 end subroutine chkin9

m_ddb/ddb_bcast [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_bcast

FUNCTION

  MPI broadcast all types for the ddb_type structure

INPUTS

   master=Rank of Master
   comm=MPI communicator

SIDE EFFECTS

   Ddb<type(ddb_type)>= Input if node is master, other nodes returns with a completely initialized instance.

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

425 subroutine ddb_bcast(Ddb, master, comm)
426 
427 
428 !This section has been created automatically by the script Abilint (TD).
429 !Do not modify the following lines by hand.
430 #undef ABI_FUNC
431 #define ABI_FUNC 'ddb_bcast'
432 !End of the abilint section
433 
434  implicit none
435 
436 !Arguments ------------------------------------
437 !array
438  type(ddb_type),intent(inout) :: ddb
439  integer, intent(in) :: master,comm
440 
441 !Local variables-------------------------------
442 !scalars
443  integer :: ierr
444 ! *************************************************************************
445 
446  if (xmpi_comm_size(comm) == 1) return
447 
448  DBG_ENTER("COLL")
449 
450  ! Transmit dimensions and static variables.
451  call xmpi_bcast(ddb%msize, master, comm, ierr)
452  call xmpi_bcast(ddb%nblok, master, comm, ierr)
453  call xmpi_bcast(ddb%natom, master, comm, ierr)
454  call xmpi_bcast(ddb%ntypat, master, comm, ierr)
455 
456  call xmpi_bcast(ddb%occopt, master, comm, ierr)
457  call xmpi_bcast(ddb%prtvol, master, comm, ierr)
458 
459  !real
460  call xmpi_bcast(ddb%rprim, master, comm, ierr)
461  call xmpi_bcast(ddb%gprim, master, comm, ierr)
462  call xmpi_bcast(ddb%acell, master, comm, ierr)
463 
464  ! Allocate arrays on the other nodes.
465  if (xmpi_comm_rank(comm) /= master) then
466    call ddb_malloc(ddb,ddb%msize,ddb%nblok,ddb%natom,ddb%ntypat)
467  end if
468 
469  call xmpi_bcast(ddb%flg, master, comm, ierr)
470  call xmpi_bcast(ddb%typ, master, comm, ierr)
471  call xmpi_bcast(ddb%amu, master, comm, ierr)
472  call xmpi_bcast(ddb%nrm, master, comm, ierr)
473  call xmpi_bcast(ddb%qpt, master, comm, ierr)
474  call xmpi_bcast(ddb%val, master, comm, ierr)
475 
476  DBG_EXIT("COLL")
477 
478 end subroutine ddb_bcast

m_ddb/ddb_copy [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_copy

FUNCTION

  Create object and copy all types for the ddb_type structure

PARENTS

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

298 subroutine ddb_copy(iddb, oddb)
299 
300 
301 !This section has been created automatically by the script Abilint (TD).
302 !Do not modify the following lines by hand.
303 #undef ABI_FUNC
304 #define ABI_FUNC 'ddb_copy'
305 !End of the abilint section
306 
307  implicit none
308 
309 !Arguments ------------------------------------
310 !array
311  type(ddb_type),intent(in) :: iddb
312  type(ddb_type),intent(out) :: oddb
313 
314 ! ************************************************************************
315 
316  ! Copy dimensions and static variables.
317  oddb%msize = iddb%msize
318  oddb%mpert = iddb%mpert
319  oddb%nblok = iddb%nblok
320  oddb%natom = iddb%natom
321  oddb%ntypat = iddb%ntypat
322  oddb%occopt = iddb%occopt
323  oddb%prtvol = iddb%prtvol
324 
325  oddb%rprim = iddb%rprim
326  oddb%gprim = iddb%gprim
327  oddb%acell = iddb%acell
328 
329  ! Allocate and copy the allocatable arrays.
330  call alloc_copy(iddb%flg,oddb%flg)
331  call alloc_copy(iddb%typ,oddb%typ)
332  call alloc_copy(iddb%amu,oddb%amu)
333  call alloc_copy(iddb%nrm,oddb%nrm)
334  call alloc_copy(iddb%qpt,oddb%qpt)
335  call alloc_copy(iddb%val,oddb%val)
336 
337 end subroutine ddb_copy

m_ddb/ddb_diagoq [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_diagoq

FUNCTION

  Compute the phonon frequencies at the specified q-point by performing
  a direct diagonalization of the dynamical matrix. The q-point **MUST** be
  one the points stored in the DDB file.

INPUTS

  ddb<type(ddb_type)>=Object storing the DDB results.
  crystal<type(crystal_t)> = Information on the crystalline structure.
  asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
  symdynmat=If equal to 1, the dynamical matrix is symmetrized in dfpt_phfrq before the diagonalization.
  rftyp  = 1 if non-stationary block
           2 if stationary block
           3 if third order derivatives
  qpt(3)=q-point in reduced coordinates.

OUTPUT

  phfrq(3*crystal%natom)=Phonon frequencies in Hartree
  displ_cart(2,3*%natom,3*%natom)=Phonon displacement in Cartesian coordinates
  [out_eigvec(2*3*natom*3*natom) = The igenvectors of the dynamical matrix.
  [out_displ_red(2*3*natom*3*natom) = The displacement in reduced coordinates.

PARENTS

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

2856 subroutine ddb_diagoq(ddb, crystal, qpt, asrq0, symdynmat, rftyp, phfrq, displ_cart, &
2857                       out_eigvec,out_displ_red)   ! Optional [out]
2858 
2859 
2860 !This section has been created automatically by the script Abilint (TD).
2861 !Do not modify the following lines by hand.
2862 #undef ABI_FUNC
2863 #define ABI_FUNC 'ddb_diagoq'
2864 !End of the abilint section
2865 
2866  implicit none
2867 
2868 !Arguments ------------------------------------
2869 !scalars
2870  integer,intent(in) :: rftyp,symdynmat
2871  type(ddb_type),intent(in) :: ddb
2872  type(asrq0_t),intent(inout) :: asrq0
2873  type(crystal_t),intent(in) :: crystal
2874 !arrays
2875  real(dp),intent(in) :: qpt(3)
2876  real(dp),intent(out) :: displ_cart(2,3,crystal%natom,3,crystal%natom)
2877  real(dp),intent(out) :: phfrq(3*crystal%natom)
2878  !real(dp),optional,intent(out) :: out_d2cart(2,3*crystal%natom,3*crystal%natom)
2879  real(dp),optional,intent(out) :: out_eigvec(2,3,crystal%natom,3*crystal%natom)
2880  real(dp),optional,intent(out) :: out_displ_red(2,3,crystal%natom,3*crystal%natom)
2881 
2882 !Local variables-------------------------------
2883  integer :: iblok,natom
2884 !arrays
2885  integer :: rfphon(4),rfelfd(4),rfstrs(4)
2886  real(dp) :: qphnrm(3), qphon_padded(3,3),d2cart(2,ddb%msize),my_qpt(3)
2887  real(dp) :: eigvec(2,3,crystal%natom,3*crystal%natom),eigval(3*crystal%natom)
2888 
2889 ! ************************************************************************
2890 
2891  ! Use my_qpt because dfpt_phfrq can change the q-point (very bad design)
2892  qphnrm = one; my_qpt = qpt
2893 
2894  ! Look for the information in the DDB (no interpolation here!)
2895  rfphon(1:2)=1
2896  rfelfd(1:2)=0
2897  rfstrs(1:2)=0
2898  qphon_padded = zero
2899  qphon_padded(:,1) = qpt
2900  natom = crystal%natom
2901 
2902  call gtblk9(ddb,iblok,qphon_padded,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2903  if (iblok == 0) then
2904    MSG_ERROR(sjoin("Cannot find q-point ", ktoa(qpt)," in DDB file"))
2905  end if
2906 
2907  ! Copy the dynamical matrix in d2cart
2908  d2cart(:,1:ddb%msize) = ddb%val(:,:,iblok)
2909 
2910  ! Eventually impose the acoustic sum rule based on previously calculated d2asr
2911  call asrq0_apply(asrq0, natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart)
2912 
2913  ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
2914  call dfpt_phfrq(ddb%amu,displ_cart,d2cart,eigval,eigvec,crystal%indsym,&
2915 &  ddb%mpert,crystal%nsym,natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),my_qpt,&
2916 &  crystal%rprimd,symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol)
2917 
2918  ! Return the dynamical matrix and the eigenvector for this q-point
2919  !if (present(out_d2cart)) out_d2cart = d2cart(:,:3*natom,:3*natom)
2920  if (present(out_eigvec)) out_eigvec = eigvec
2921 
2922  ! Return phonon displacement in reduced coordinates.
2923  if (present(out_displ_red)) call phdispl_cart2red(natom, crystal%gprimd, displ_cart, out_displ_red)
2924 
2925 end subroutine ddb_diagoq

m_ddb/ddb_free [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_free

FUNCTION

  Clean and deallocate types for the ddb_type structure

PARENTS

      anaddb,ddb_interpolate,dfpt_looppert,dfptnl_doutput,eph,gstate
      m_effective_potential_file,m_gruneisen,mblktyp1,mblktyp5,thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

241 subroutine ddb_free(ddb)
242 
243 
244 !This section has been created automatically by the script Abilint (TD).
245 !Do not modify the following lines by hand.
246 #undef ABI_FUNC
247 #define ABI_FUNC 'ddb_free'
248 !End of the abilint section
249 
250  implicit none
251 
252 !Arguments ------------------------------------
253  type(ddb_type),intent(inout) :: ddb
254 
255 ! ************************************************************************
256 
257  !integer
258  if (allocated(ddb%flg))  then
259    ABI_FREE(ddb%flg)
260  end if
261  if (allocated(ddb%typ))  then
262    ABI_FREE(ddb%typ)
263  end if
264 
265  ! real
266  if (allocated(ddb%amu))  then
267    ABI_FREE(ddb%amu)
268  end if
269  if (allocated(ddb%nrm))  then
270    ABI_FREE(ddb%nrm)
271  end if
272  if (allocated(ddb%qpt))  then
273    ABI_FREE(ddb%qpt)
274  end if
275  if (allocated(ddb%val))  then
276    ABI_FREE(ddb%val)
277  end if
278 
279 end subroutine ddb_free

m_ddb/ddb_from_file [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_from_file

FUNCTION

  This subroutine reads data from the DDB file and constructs an instance of ddb_type
  It also returns an instance of crystal_t with the crystalline structure reported in the DDB file

INPUTS

  filename=DDB filename.
  brav = 1 or -1 -> simple lattice; 2 -> face-centered cubic;
         3 -> body-centered lattice; 4 -> hexagonal lattice (D6h)
  natom=Number of atoms in the unit cell
  atifc(natom)=list of the atom ifc to be analysed
  natifc = number of atoms for which the analysis of ifc is done
  comm=MPI communicator.
  [prtvol] = Verbosity level

OUTPUT

  ddb<type(ddb_type)>=Object storing the DDB results.
  Crystal<type(crystal_t)>=Crystal structure parameters
  atifc(natom) =  atifc(ia) equals 1 if the analysis of ifc
    has to be done for atom ia; otherwise 0.

TODO

   Sorry for the presence of natom, natifc and atifc.
   They are needed for legacy code!

PARENTS

      anaddb,dfpt_looppert,eph,m_effective_potential_file,m_gruneisen

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1633 subroutine ddb_from_file(ddb,filename,brav,natom,natifc,atifc,crystal,comm,prtvol)
1634 
1635 
1636 !This section has been created automatically by the script Abilint (TD).
1637 !Do not modify the following lines by hand.
1638 #undef ABI_FUNC
1639 #define ABI_FUNC 'ddb_from_file'
1640 !End of the abilint section
1641 
1642  implicit none
1643 
1644 !Arguments ------------------------------------
1645 !scalars
1646  integer,intent(in) :: comm,brav,natom,natifc
1647  integer,optional,intent(in) :: prtvol
1648  character(len=*),intent(in) :: filename
1649  type(crystal_t),intent(out) :: Crystal
1650  type(ddb_type),intent(inout) :: ddb
1651 !array
1652  integer,intent(inout) :: atifc(natom)
1653 
1654 !Local variables-------------------------------
1655 !scalars
1656  integer,parameter :: master=0
1657  integer :: ierr,ii,msym,dimekb,lmnmax,mband,nkpt,ntypat,nsym,usepaw
1658  integer :: mtyp,mpert,msize,ddb_natom,nblok,occopt,timrev,space_group,npsp,ddbun
1659  real(dp) :: factor,ucvol
1660  logical :: use_antiferro
1661  type(ddb_hdr_type) :: ddb_hdr
1662 !arrays
1663  integer,allocatable :: symrec(:,:,:),symrel(:,:,:),symafm(:),indsym(:,:,:),typat(:)
1664  real(dp) :: acell(3),gmet(3,3),gprim(3,3),rmet(3,3),rprim(3,3),rprimd(3,3)
1665  real(dp),allocatable :: amu(:),xcart(:),xred(:,:),zion(:),znucl(:),tnons(:,:)
1666  character(len=132),allocatable :: title(:)
1667  character(len=500) :: message
1668 
1669 ! ************************************************************************
1670 
1671  DBG_ENTER("COLL")
1672 
1673 ! Must read natom from the DDB before being able to allocate some arrays needed for invars9
1674  ddbun = get_unit()
1675  call ddb_hdr_open_read(ddb_hdr,filename,ddbun,DDB_VERSION, comm=comm, dimonly=1)
1676 
1677  nblok = ddb_hdr%nblok
1678  mtyp = ddb_hdr%mblktyp
1679  msym = ddb_hdr%msym
1680  ddb_natom = ddb_hdr%natom
1681  ntypat = ddb_hdr%ntypat
1682  mband = ddb_hdr%mband
1683  nkpt = ddb_hdr%nkpt
1684  usepaw = ddb_hdr%usepaw
1685  dimekb = ddb_hdr%psps%dimekb
1686  lmnmax = ddb_hdr%psps%lmnmax
1687 
1688  ! JWZ occopt was used below before being initialized 13 April 2018
1689  occopt = ddb_hdr%occopt
1690  call ddb_hdr_free(ddb_hdr)
1691 
1692  if (ddb_natom /= natom) then
1693    MSG_ERROR(sjoin("input natom:",itoa(natom),"does not agree with DDB value:",itoa(natom)))
1694  end if
1695 
1696  mpert=natom+6
1697  msize=3*mpert*3*mpert; if (mtyp==3) msize=msize*3*mpert
1698 
1699  ! Allocate arrays depending on msym (which is actually fixed to nsym inside inprep8)
1700  ABI_MALLOC(symrel,(3,3,msym))
1701  ABI_MALLOC(symafm,(msym))
1702  ABI_MALLOC(tnons,(3,msym))
1703 
1704  ABI_MALLOC(typat,(natom))
1705  ABI_MALLOC(xred,(3,natom))
1706  ABI_MALLOC(zion,(ntypat))
1707  ABI_MALLOC(znucl,(ntypat))
1708 
1709  ! Master reads and then broadcasts data.
1710  if (xmpi_comm_rank(comm) == master) then
1711 
1712    ABI_MALLOC(symrec,(3,3,msym))
1713    ABI_MALLOC(indsym,(4,msym,natom))
1714    ABI_MALLOC(xcart,(3*natom))
1715    ABI_MALLOC(amu,(ntypat))
1716 
1717    ddbun = get_unit() ! FIXME: The treatment of the unit number in rdddb9 is ugly!
1718 
1719    call ddb_malloc(ddb,msize,nblok,natom,ntypat)
1720 
1721    call rdddb9(acell,atifc,amu,ddb,&
1722 &   ddbun,filename,gmet,gprim,indsym,ab_out,&
1723 &   mband,mpert,msize,msym,&
1724 &   natifc,ddb_natom,nkpt,nsym,ntypat,&
1725 &   rmet,rprim,symrec,symrel,symafm,&
1726 &   tnons,typat,ucvol,xcart,xred,zion,znucl)
1727 
1728    close(ddbun)
1729 
1730    ABI_FREE(symrec)
1731    ABI_FREE(indsym)
1732    ABI_FREE(xcart)
1733 
1734    ! Renormalize rprim to possibly satisfy the constraint abs(rprim(1,2))=half when abs(brav)/=1
1735    ! This section is needed to preserver the behaviour of the old implementation.
1736    if (abs(brav)/=1 .and. abs(abs(rprim(1,2))-half)>tol10) then
1737      if(abs(rprim(1,2))<tol6)then
1738        write(message, '(a,i0,7a)' )&
1739 &       'The input DDB value of brav is ',brav,',',ch10,&
1740 &       'and the one of rprim(1,2) is zero.',ch10,&
1741 &       'These are incompatible',ch10,&
1742 &       'Action: check the value of brav and rprim(1,2) in your DDB.'
1743        MSG_ERROR(message)
1744      end if
1745      factor=abs(rprim(1,2))*two
1746      acell(:)=acell(:)*factor
1747      rprim(:,:)=rprim(:,:)/factor
1748      gprim(:,:)=gprim(:,:)*factor
1749    end if
1750 
1751    ! Save variables needed to call legacy code.
1752    ddb%acell = acell
1753    ddb%rprim = rprim
1754    ddb%gprim = gprim
1755 
1756    ! Other useful quantities.
1757    ! 2 is to preserve the old behaviour
1758    ddb%prtvol = 2; if (present(prtvol)) ddb%prtvol = prtvol
1759    ddb%occopt = occopt
1760    ddb%amu = amu
1761    ABI_FREE(amu)
1762 
1763    ! Now the whole DDB is in central memory, contained in the array ddb%val(2,msize,nblok).
1764    ! The data is contained in the four arrays
1765    !   ddb%flg(msize,nblok) : blok flag for each element
1766    !   ddb%qpt(9,nblok)     : blok wavevector (unnormalized)
1767    !   ddb%nrm(3,nblok)     : blok wavevector normalization
1768    !   ddb%typ(nblok)       : blok type
1769  end if
1770 
1771  if (xmpi_comm_size(comm) > 1) then
1772    call ddb_bcast (ddb, master, comm)
1773    call xmpi_bcast(atifc, master, comm, ierr)
1774    call xmpi_bcast(nsym, master, comm, ierr)
1775    call xmpi_bcast(symrel, master, comm, ierr)
1776    call xmpi_bcast(symafm, master, comm, ierr)
1777    call xmpi_bcast(typat, master, comm, ierr)
1778    call xmpi_bcast(acell, master, comm, ierr)
1779    call xmpi_bcast(occopt, master, comm, ierr)
1780    call xmpi_bcast(gprim, master, comm, ierr)
1781    call xmpi_bcast(rprim, master, comm, ierr)
1782    call xmpi_bcast(tnons, master, comm, ierr)
1783    call xmpi_bcast(xred, master, comm, ierr)
1784    call xmpi_bcast(zion, master, comm, ierr)
1785    call xmpi_bcast(znucl, master, comm, ierr)
1786  end if
1787 
1788 !Initialize crystal_t object.
1789  call mkrdim(acell,rprim,rprimd)
1790 
1791 !FIXME: These variables are hardcoded
1792  npsp = ntypat; space_group = 0; timrev = 2
1793  use_antiferro=.FALSE. !;  use_antiferro=(nspden==2.and.nsppol==1)
1794  ABI_MALLOC(title, (ntypat))
1795 
1796  do ii=1,ntypat
1797    write(title(ii),'(a,i0)')"No title for typat ",ii
1798  end do
1799 
1800 !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
1801  call crystal_init(ddb%amu,Crystal,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,&
1802 &  zion,znucl,timrev,use_antiferro,.FALSE.,title,&
1803 &  symrel=symrel(:,:,1:nsym),tnons=tnons(:,1:nsym),symafm=symafm(1:nsym))
1804 
1805  ABI_FREE(title)
1806  ABI_FREE(symrel)
1807  ABI_FREE(symafm)
1808  ABI_FREE(tnons)
1809  ABI_FREE(typat)
1810  ABI_FREE(xred)
1811  ABI_FREE(zion)
1812  ABI_FREE(znucl)
1813 
1814  DBG_EXIT("COLL")
1815 
1816 end subroutine ddb_from_file

m_ddb/ddb_get_asrq0 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_get_asrq0

FUNCTION

  In case the interatomic forces are not calculated, the
  ASR-correction has to be determined here from the Dynamical matrix at Gamma.
  In case the DDB does not contain this information, the subroutine returns iblok=0
  %d2asr is initialized and set to zero to preserve the old behaviour.

INPUTS

  asr=Input variable selecting the method for the ASR
  rftyp  = 1 if non-stationary block
           2 if stationary block
           3 if third order derivatives
  xcart(3,ddb%atom)=Cartesian coordinates of the atoms.

SIDE EFFECTS

  ddb<type(ddb_type)>= Database with the derivates. The routine does not change it
  except when asr is in [3,4]. TODO This should not happen.

OUTPUT

 asrq0<asrq0_t>
   iblok= is set to 0 if the Gamma block is not found

PARENTS

CHILDREN

SOURCE

2732 type(asrq0_t) function ddb_get_asrq0(ddb, asr, rftyp, xcart) result(asrq0)
2733 
2734 
2735 !This section has been created automatically by the script Abilint (TD).
2736 !Do not modify the following lines by hand.
2737 #undef ABI_FUNC
2738 #define ABI_FUNC 'ddb_get_asrq0'
2739 !End of the abilint section
2740 
2741  implicit none
2742 
2743 !Arguments ------------------------------------
2744 !scalars
2745  integer,intent(in) :: asr,rftyp
2746  type(ddb_type),intent(inout) :: ddb
2747 !arrays
2748  real(dp),intent(in) :: xcart(3,ddb%natom)
2749 
2750 !Local variables-------------------------------
2751 !scalars
2752  integer :: dims,iblok
2753  !character(len=500) :: msg
2754 !arrays
2755  integer :: rfelfd(4),rfphon(4),rfstrs(4)
2756  real(dp) :: qphnrm(3),qphon(3,3)
2757  real(dp),allocatable :: d2asr_res(:,:,:,:,:),d2cart(:,:)
2758 
2759 ! ************************************************************************
2760 
2761  asrq0%asr = asr; asrq0%natom = ddb%natom
2762 
2763  ! Find the Gamma block in the DDB (no need for E-field entries)
2764  qphon(:,1)=zero
2765  qphnrm(1)=zero
2766  rfphon(1:2)=1
2767  rfelfd(:)=0
2768  rfstrs(:)=0
2769 
2770  call gtblk9(ddb,asrq0%iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2771  ! this is to maintain the old behaviour in which the arrays where allocated and set to zero in anaddb.
2772  ABI_MALLOC(asrq0%d2asr, (2,3,ddb%natom,3,ddb%natom))
2773  asrq0%d2asr = zero
2774 
2775  if (asrq0%iblok == 0) return
2776  iblok = asrq0%iblok
2777 
2778  select case (asrq0%asr)
2779  case (0)
2780    continue
2781 
2782  case (1,2)
2783    call asria_calc(asr,asrq0%d2asr,ddb%val(:,:,iblok),ddb%mpert,ddb%natom)
2784 
2785  case (3,4)
2786    ! Rotational invariance for 1D and 0D systems
2787    ! Compute uinvers, vtinvers and singular matrices.
2788    dims = 3*ddb%natom*(3*ddb%natom-1) / 2
2789    ABI_CALLOC(asrq0%uinvers, (dims, dims))
2790    ABI_CALLOC(asrq0%vtinvers,(dims, dims))
2791    ABI_CALLOC(asrq0%singular, (dims))
2792 
2793    call asrprs(asr,1,3,asrq0%uinvers,asrq0%vtinvers,asrq0%singular,&
2794      ddb%val(:,:,iblok),ddb%mpert,ddb%natom,xcart)
2795 
2796  case (5)
2797    ! d2cart is a temp variable here
2798    ABI_MALLOC(d2cart,(2,ddb%msize))
2799    d2cart = ddb%val(:,:,iblok)
2800    ! calculate diagonal correction
2801    call asria_calc(2,asrq0%d2asr,d2cart,ddb%mpert,ddb%natom)
2802    ! apply diagonal correction
2803    call asria_corr(2,asrq0%d2asr,d2cart,ddb%mpert,ddb%natom)
2804    ! hermitianize
2805    call mkherm(d2cart,3*ddb%mpert)
2806    ! remove remaining ASR rupture due to Hermitianization
2807    ABI_MALLOC(d2asr_res,(2,3,ddb%natom,3,ddb%natom))
2808    call asria_calc(asr,d2asr_res,d2cart,ddb%mpert,ddb%natom)
2809    ! full correction is sum of both
2810    asrq0%d2asr = asrq0%d2asr + d2asr_res
2811 
2812    ABI_FREE(d2cart)
2813    ABI_FREE(d2asr_res)
2814 
2815  case default
2816    MSG_ERROR(sjoin("Wrong value for asr:", itoa(asr)))
2817  end select
2818 
2819 end function ddb_get_asrq0

m_ddb/ddb_get_dchidet [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_get_dchidet

FUNCTION

 Reads the non-linear optical susceptibility tensor and the
 first-order change in the linear dielectric susceptibility

INPUTS

 ddb<type(ddb_type)>=Derivative Database.
 ramansr= if /= 0, impose sum rule on first-order derivatives
                   of the electronic susceptibility with respect
                   to atomic displacements
 nlflag= if =3, only the non-linear optical susceptibilities is computed

OUTPUT

 dchide(3,3,3) = non-linear optical coefficients
 dchidt(natom,3,3,3) = first-order change of the electronic dielectric
   tensor induced by an individual atomic displacement
 iblok=Index of the block containing the data. 0 if block is not found.
   The caller should check the returned value.

PARENTS

CHILDREN

SOURCE

2647 integer function ddb_get_dchidet(ddb,ramansr,nlflag,dchide,dchidt) result(iblok)
2648 
2649 
2650 !This section has been created automatically by the script Abilint (TD).
2651 !Do not modify the following lines by hand.
2652 #undef ABI_FUNC
2653 #define ABI_FUNC 'ddb_get_dchidet'
2654 !End of the abilint section
2655 
2656  implicit none
2657 
2658 !Arguments -------------------------------
2659 !scalars
2660  integer,intent(in) :: ramansr, nlflag
2661  type(ddb_type),intent(in) :: ddb
2662 !arrays
2663  real(dp),intent(out) :: dchide(3,3,3),dchidt(ddb%natom,3,3,3)
2664 
2665 !Local variables -------------------------
2666 !scalars
2667  integer :: rftyp
2668 !arrays
2669  integer :: rfelfd(4),rfphon(4),rfstrs(4)
2670  real(dp) :: qphnrm(3),qphon(3,3)
2671 
2672 ! *********************************************************************
2673 
2674  qphon(:,:) = zero
2675  qphnrm(:)  = one
2676 ! rfphon(1)  = 1 ; rfphon(2:3) = 0
2677  rfelfd(:)  = 2
2678  rfstrs(:)  = 0
2679  rftyp = 3
2680 
2681  if (nlflag < 3) then
2682    rfphon(1)  = 1 ; rfphon(2:3) = 0
2683  else
2684    rfphon(1)  = 0 ; rfphon(2:3) = 0
2685  end if
2686 
2687  call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2688 
2689  if (iblok /= 0) then
2690    call dtchi(ddb%val(:,:,iblok),dchide,dchidt,ddb%mpert,ddb%natom,ramansr,nlflag)
2691  else
2692    ! Let the caller handle the error.
2693    dchide = huge(one); dchidt = huge(one)
2694  end if
2695 
2696 end function ddb_get_dchidet

m_ddb/ddb_get_dielt [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_get_dielt

FUNCTION

 Reads the Dielectric Tensor from the DDB file

INPUTS

  ddb<type(ddb_type)>=Derivative database.
  rftyp  = 1 if non-stationary block
           2 if stationary block
           3 if third order derivatives

OUTPUT

  dielt(3,3) = Macroscopic dielectric tensor
  iblok=Index of the block containing the data. 0 if block is not found.

NOTES

  dielt is initialized to one_3D if the derivatives are not available in the DDB file.

PARENTS

CHILDREN

SOURCE

2552 integer function ddb_get_dielt(ddb,rftyp,dielt) result(iblok)
2553 
2554 
2555 !This section has been created automatically by the script Abilint (TD).
2556 !Do not modify the following lines by hand.
2557 #undef ABI_FUNC
2558 #define ABI_FUNC 'ddb_get_dielt'
2559 !End of the abilint section
2560 
2561  implicit none
2562 
2563 !Arguments -------------------------------
2564 !scalars
2565  integer,intent(in) :: rftyp
2566  type(ddb_type),intent(in) :: ddb
2567 !arrays
2568  real(dp),intent(out) :: dielt(3,3)
2569 
2570 !Local variables -------------------------
2571 !scalars
2572  integer :: mpert
2573  character(len=1000) :: message
2574 !arrays
2575  integer :: rfelfd(4),rfphon(4),rfstrs(4)
2576  real(dp) :: qphnrm(3),qphon(3,3)
2577  real(dp),allocatable :: tmpval(:,:,:,:)
2578 
2579 ! *********************************************************************
2580 
2581  ! Look for the Gamma Block in the DDB
2582  qphon(:,1)=zero
2583  qphnrm(1)=zero
2584  rfphon(1:2)=0
2585  rfelfd(1:2)=2
2586  rfstrs(1:2)=0
2587 
2588  call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2589 
2590  ! Read the dielectric tensor only if the Gamma-block was found in the DDB
2591  ! In case it was not found, iblok = 0
2592  dielt=zero; dielt(1,1)=one; dielt(2,2)=one; dielt(3,3)=one
2593 
2594  if (iblok/=0) then
2595    !Extration of dielectric tensor
2596    mpert = ddb%mpert
2597 
2598    ABI_MALLOC(tmpval,(3,mpert,3,mpert))
2599    tmpval(:,:,:,:) = reshape(ddb%val(1,:,iblok), shape = (/3,mpert,3,mpert/))
2600    dielt=tmpval(1:3,ddb%natom+2,1:3,ddb%natom+2)
2601 
2602    write(message,'(a,3es16.6,3es16.6,3es16.6)' )&
2603 &   ' Dielectric Tensor ',&
2604 &   dielt(1,1),dielt(1,2),dielt(1,3),&
2605 &   dielt(2,1),dielt(2,2),dielt(2,3),&
2606 &   dielt(3,1),dielt(3,2),dielt(3,3)
2607 
2608    call wrtout(std_out,message,'COLL')
2609 
2610    ABI_FREE(tmpval)
2611  end if ! iblok not found
2612 
2613 end function ddb_get_dielt

m_ddb/ddb_get_dielt_zeff [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_get_dielt_zeff

FUNCTION

 Reads the Dielectric Tensor and the Effective Charges from the DDB file
 Impose the charge neutrality on the effective charges and eventually select some parts of the effective charges

INPUTS

  Crystal<type(crystal_t)>=Crystal structure parameters
  rftyp  = 1 if non-stationary block
           2 if stationary block
           3 if third order derivatives
  chneut=(0 => no ASR, 1 => equal repartition,2 => weighted repartition )
  selectz=selection of some parts of the effective charge tensor attached to one atom.
    (0=> no selection, 1=> trace only, 2=> symmetric part only)                       !!

SIDE EFFECTS

  ddb<type(ddb_type)>=
    The block with the effective charges is modified if charge neutrality is imposed.

OUTPUT

  dielt(3,3) = Macroscopic dielectric tensor
  zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
  iblok=Index of the block containing the data. 0 if block is not found.

NOTES

  dielt and zeff are initialized to one_3D and zero if the derivatives are not available in the DDB file.

PARENTS

CHILDREN

SOURCE

2456 integer function ddb_get_dielt_zeff(ddb,crystal,rftyp,chneut,selectz,dielt,zeff) result(iblok)
2457 
2458 
2459 !This section has been created automatically by the script Abilint (TD).
2460 !Do not modify the following lines by hand.
2461 #undef ABI_FUNC
2462 #define ABI_FUNC 'ddb_get_dielt_zeff'
2463 !End of the abilint section
2464 
2465  implicit none
2466 
2467 !Arguments -------------------------------
2468 !scalars
2469  integer,intent(in) :: rftyp,chneut,selectz
2470  type(ddb_type),intent(inout) :: ddb
2471  type(crystal_t),intent(in) :: crystal
2472 !arrays
2473  real(dp),intent(out) :: dielt(3,3),zeff(3,3,crystal%natom)
2474 
2475 !Local variables -------------------------
2476 !scalars
2477  integer :: ii
2478  character(len=500) :: message
2479 !arrays
2480  integer :: rfelfd(4),rfphon(4),rfstrs(4)
2481  real(dp) :: qphnrm(3),qphon(3,3)
2482 
2483 ! *********************************************************************
2484 
2485  ! Look for the Gamma Block in the DDB
2486  qphon(:,1)=zero
2487  qphnrm(1)=zero
2488  rfphon(1:2)=1
2489  rfelfd(1:2)=2
2490  rfstrs(1:2)=0
2491 
2492  !write(std_out,*)"ddb%mpert",ddb%mpert
2493 
2494  call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2495 
2496  ! Compute effective charges and dielectric tensor only if the Gamma-blok was found in the DDB
2497  ! In case it was not found, iblok = 0
2498  zeff=zero; dielt=zero; dielt(1,1)=one; dielt(2,2)=one; dielt(3,3)=one
2499 
2500  if (iblok/=0) then
2501    write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,&
2502 &   ' Dielectric Tensor and Effective Charges ',ch10
2503    call wrtout(std_out,message,'COLL')
2504    call wrtout(ab_out,message,'COLL')
2505 
2506    ! Make the imaginary part of the Gamma block vanish
2507    write(message, '(a,a,a,a,a)'  ) ch10,&
2508 &   ' anaddb : Zero the imaginary part of the Dynamical Matrix at Gamma,',ch10,&
2509 &   '   and impose the ASR on the effective charges ',ch10
2510    call wrtout(std_out,message,'COLL')
2511    call wrtout(ab_out,message,'COLL')
2512 
2513    ! Impose the charge neutrality on the effective charges and eventually select some parts of the effective charges
2514    call chneu9(chneut,ddb%val(:,:,iblok),ddb%mpert,ddb%natom,ddb%ntypat,selectz,Crystal%typat,Crystal%zion)
2515 
2516    ! Extraction of the dielectric tensor and the effective charges
2517    call dtech9(ddb%val,dielt,iblok,ddb%mpert,ddb%natom,ddb%nblok,zeff)
2518  end if ! iblok not found
2519 
2520 end function ddb_get_dielt_zeff

m_ddb/ddb_get_etotal [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

  ddb_get_etotal

FUNCTION

  Read the GS total energy from the DDB file.

INPUTS

  ddb<type(ddb_type)>=Derivative database.

OUTPUT

  etotal=GS Total energy in Hartree
  iblock=Index of the block in the DDB file. 0 if not found.

PARENTS

CHILDREN

SOURCE

2374 integer function ddb_get_etotal(ddb,etotal) result(iblok)
2375 
2376 
2377 !This section has been created automatically by the script Abilint (TD).
2378 !Do not modify the following lines by hand.
2379 #undef ABI_FUNC
2380 #define ABI_FUNC 'ddb_get_etotal'
2381 !End of the abilint section
2382 
2383  implicit none
2384 
2385 !Arguments -------------------------------
2386 !scalars
2387  real(dp),intent(out) :: etotal
2388  type(ddb_type),intent(in) :: ddb
2389 
2390 !Local variables -------------------------
2391 !scalars
2392  integer :: rftyp
2393 !arrays
2394  integer :: rfelfd(4),rfphon(4),rfstrs(4)
2395  real(dp) :: qphnrm(3),qphon(3,3)
2396 
2397 ! *********************************************************************
2398 
2399  ! Extract the block with the total energy
2400  qphon(:,:) = zero
2401  qphnrm(:) = zero
2402  rfphon(:) = 0
2403  rfelfd(:) = 0
2404  rfstrs(:) = 0
2405  rftyp = 0
2406 
2407  call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2408 
2409  if (iblok /= 0) then
2410    etotal = ddb%val(1,1,iblok)
2411  else
2412    etotal = huge(one)
2413  end if
2414 
2415 end function ddb_get_etotal

m_ddb/ddb_malloc [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_malloc

FUNCTION

INPUTS

OUTPUT

PARENTS

      ddb_interpolate,dfptnl_doutput,gstate,m_ddb,mblktyp1,mblktyp5,thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

360 subroutine ddb_malloc(ddb,msize,nblok,natom,ntypat)
361 
362 
363 !This section has been created automatically by the script Abilint (TD).
364 !Do not modify the following lines by hand.
365 #undef ABI_FUNC
366 #define ABI_FUNC 'ddb_malloc'
367 !End of the abilint section
368 
369  implicit none
370 
371 !Arguments ------------------------------------
372 !array
373  integer,intent(in) :: msize,nblok,natom,ntypat
374  type(ddb_type),intent(inout) :: ddb
375 
376 ! ************************************************************************
377 
378  ! FIXME
379  ! This is done in rdddb9 but only by the master node!
380  ! Should rationalize the different steps
381  ddb%msize = msize
382  ddb%nblok = nblok
383  ddb%natom = natom
384  ddb%mpert = natom+6
385  ddb%ntypat = ntypat
386 
387  ! integer
388  ABI_CALLOC(ddb%flg,(msize,nblok))
389  ABI_CALLOC(ddb%typ,(nblok))
390 
391  ! real
392  ABI_MALLOC(ddb%amu,(ntypat))
393  ABI_MALLOC(ddb%nrm,(3,nblok))
394  ABI_MALLOC(ddb%qpt,(9,nblok))
395  ABI_MALLOC(ddb%val,(2,msize,nblok))
396  ddb%val = huge(one)
397 
398 end subroutine ddb_malloc

m_ddb/ddb_to_dtset [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_to_dtset

FUNCTION

   Initialize a dataset object from ddb.

INPUTS

OUTPUT

PARENTS

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

3393 subroutine ddb_to_dtset(comm,dtset,filename,psps)
3394 
3395 
3396 !This section has been created automatically by the script Abilint (TD).
3397 !Do not modify the following lines by hand.
3398 #undef ABI_FUNC
3399 #define ABI_FUNC 'ddb_to_dtset'
3400 !End of the abilint section
3401 
3402  implicit none
3403 
3404  !Arguments ------------------------------------
3405  integer,intent(in) :: comm
3406  type(dataset_type),intent(inout) :: dtset
3407  type(pseudopotential_type),intent(inout) :: psps
3408  ! type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
3409  character(len=*),intent(in) :: filename
3410  !Local variables -------------------------
3411  integer :: mxnimage,ddbun
3412 !integer :: ii, nn
3413  type(ddb_hdr_type) :: ddb_hdr
3414 
3415 ! ************************************************************************
3416 
3417  ABI_UNUSED(psps%usepaw)
3418 
3419 !Set variables
3420  mxnimage = 1 ! Only 1 image in the DDB
3421 
3422 ! Must read natom from the DDB before being able to allocate some arrays needed for invars9
3423  ddbun = get_unit()
3424  call ddb_hdr_open_read(ddb_hdr,filename,ddbun,DDB_VERSION,comm=comm)
3425 !close ddb file, just want to read the headers
3426  close(ddbun)
3427  dtset%ngfft = ddb_hdr%ngfft
3428 
3429 ! call psps_copy(psps, ddb_hdr%psps)
3430 
3431 ! Copy scalars from ddb
3432  dtset%natom = ddb_hdr%natom
3433  dtset%mband = ddb_hdr%mband
3434  dtset%nkpt = ddb_hdr%nkpt
3435  dtset%nsym = ddb_hdr%msym
3436  dtset%ntypat = ddb_hdr%ntypat
3437  dtset%nspden = ddb_hdr%nspden
3438  dtset%nspinor = ddb_hdr%nspinor
3439  dtset%nsppol = ddb_hdr%nsppol
3440  dtset%occopt = ddb_hdr%occopt
3441  dtset%usepaw = ddb_hdr%usepaw
3442  dtset%intxc = ddb_hdr%intxc
3443  dtset%ixc = ddb_hdr%ixc
3444  dtset%iscf = ddb_hdr%iscf
3445  dtset%dilatmx = ddb_hdr%dilatmx
3446  dtset%ecut = ddb_hdr%ecut
3447  dtset%ecutsm = ddb_hdr%ecutsm
3448  dtset%pawecutdg = ddb_hdr%pawecutdg
3449  dtset%kptnrm = ddb_hdr%kptnrm
3450  dtset%dfpt_sciss = ddb_hdr%dfpt_sciss
3451  dtset%tolwfr = 1.0_dp  ! dummy
3452  dtset%tphysel = ddb_hdr%tphysel
3453  dtset%tsmear = ddb_hdr%tsmear
3454 
3455  ! Copy arrays from ddb
3456  if (allocated(dtset%acell_orig)) then
3457    ABI_DEALLOCATE(dtset%acell_orig)
3458  end if
3459  ABI_ALLOCATE(dtset%acell_orig,(3,mxnimage))
3460  dtset%acell_orig(1:3,1) = ddb_hdr%acell(:)
3461 
3462  if (allocated(dtset%rprim_orig)) then
3463    ABI_DEALLOCATE(dtset%rprim_orig)
3464  end if
3465  ABI_ALLOCATE(dtset%rprim_orig,(3,3,mxnimage))
3466  dtset%rprim_orig(1:3,1:3,1) = ddb_hdr%rprim(:,:)
3467 
3468  if (allocated(dtset%rprimd_orig)) then
3469    ABI_DEALLOCATE(dtset%rprimd_orig)
3470  end if
3471  ABI_ALLOCATE(dtset%rprimd_orig,(3,3,mxnimage))
3472  dtset%rprimd_orig(:,1,1) = ddb_hdr%rprim(:,1) * dtset%acell_orig(1,1)
3473  dtset%rprimd_orig(:,2,1) = ddb_hdr%rprim(:,2) * dtset%acell_orig(2,1)
3474  dtset%rprimd_orig(:,3,1) = ddb_hdr%rprim(:,3) * dtset%acell_orig(3,1)
3475 
3476  if (allocated(dtset%amu_orig)) then
3477    ABI_DEALLOCATE(dtset%amu_orig)
3478  end if
3479  ABI_ALLOCATE(dtset%amu_orig,(dtset%ntypat,mxnimage))
3480  dtset%amu_orig(:,1) = ddb_hdr%amu(:)
3481 
3482  if (allocated(dtset%typat)) then
3483    ABI_DEALLOCATE(dtset%typat)
3484  end if
3485  ABI_ALLOCATE(dtset%typat,(dtset%natom))
3486  dtset%typat(:) = ddb_hdr%typat(1:ddb_hdr%matom)
3487 
3488  if (allocated(dtset%spinat)) then
3489    ABI_DEALLOCATE(dtset%spinat)
3490  end if
3491  ABI_ALLOCATE(dtset%spinat,(3,dtset%natom))
3492  dtset%spinat(:,:) = ddb_hdr%spinat(1:3,1:ddb_hdr%matom)
3493 
3494  if (allocated(dtset%xred_orig)) then
3495    ABI_DEALLOCATE(dtset%xred_orig)
3496  end if
3497  ABI_ALLOCATE(dtset%xred_orig,(3,dtset%natom,mxnimage))
3498  dtset%xred_orig(:,:,1) = ddb_hdr%xred(1:3,1:ddb_hdr%matom)
3499 
3500  if (allocated(dtset%ziontypat)) then
3501    ABI_DEALLOCATE(dtset%ziontypat)
3502  end if
3503  ABI_ALLOCATE(dtset%ziontypat,(dtset%ntypat))
3504  dtset%ziontypat(1:ddb_hdr%mtypat) = ddb_hdr%zion(1:ddb_hdr%mtypat)
3505 
3506  if (allocated(dtset%znucl)) then
3507    ABI_DEALLOCATE(dtset%znucl)
3508  end if
3509  ABI_ALLOCATE(dtset%znucl,(dtset%ntypat))
3510  dtset%znucl(:) = ddb_hdr%znucl(1:ddb_hdr%mtypat)
3511 
3512  if (allocated(dtset%nband)) then
3513    ABI_DEALLOCATE(dtset%nband)
3514  end if
3515  ABI_ALLOCATE(dtset%nband,(dtset%nkpt))
3516  dtset%nband(:) = ddb_hdr%nband(1:ddb_hdr%mkpt*ddb_hdr%nsppol)
3517 
3518  if (allocated(dtset%symafm)) then
3519    ABI_DEALLOCATE(dtset%symafm)
3520  end if
3521  ABI_ALLOCATE(dtset%symafm,(dtset%nsym))
3522  dtset%symafm(:) = ddb_hdr%symafm(1:ddb_hdr%msym)
3523 
3524  if (allocated(dtset%symrel)) then
3525    ABI_DEALLOCATE(dtset%symrel)
3526  end if
3527  ABI_ALLOCATE(dtset%symrel,(3,3,dtset%nsym))
3528  dtset%symrel(:,:,:) = ddb_hdr%symrel(1:3,1:3,1:ddb_hdr%msym)
3529 
3530  if (allocated(dtset%tnons)) then
3531    ABI_DEALLOCATE(dtset%tnons)
3532  end if
3533  ABI_ALLOCATE(dtset%tnons,(3,dtset%nsym))
3534  dtset%tnons(:,:) = ddb_hdr%tnons(1:3,1:ddb_hdr%msym)
3535 
3536  if (allocated(dtset%kpt)) then
3537    ABI_DEALLOCATE(dtset%kpt)
3538  end if
3539  ABI_ALLOCATE(dtset%kpt,(3,dtset%nkpt))
3540  dtset%kpt(:,:) = ddb_hdr%kpt(1:3,1:ddb_hdr%mkpt)
3541 
3542  if (allocated(dtset%wtk)) then
3543    ABI_DEALLOCATE(dtset%wtk)
3544  end if
3545  ABI_ALLOCATE(dtset%wtk,(dtset%nkpt))
3546  dtset%wtk(:) = ddb_hdr%wtk(1:ddb_hdr%mkpt)
3547 
3548  ! GA: I had way too much problems implementing pawtab_copy.
3549  !     The script check-libpaw would report all sorts of errors.
3550  !     Therefore, I do a cheap copy here, copying only the relevant info.
3551  !call pawtab_copy(pawtab, ddb_hdr%pawtab)
3552  ! nn=size(pawtab)
3553  ! if (nn.gt.0) then
3554  !   do ii=1,nn
3555  !     pawtab(ii)%basis_size =ddb_hdr%pawtab(ii)%basis_size
3556  !     pawtab(ii)%lmn_size =ddb_hdr%pawtab(ii)%lmn_size
3557  !     pawtab(ii)%lmn2_size =ddb_hdr%pawtab(ii)%lmn2_size
3558  !     pawtab(ii)%rpaw =ddb_hdr%pawtab(ii)%rpaw
3559  !     pawtab(ii)%rshp =ddb_hdr%pawtab(ii)%rshp
3560  !     pawtab(ii)%shape_type =ddb_hdr%pawtab(ii)%shape_type
3561  !    if (allocated(pawtab(ii)%dij0)) then
3562  !      call alloc_copy(ddb_hdr%pawtab(ii)%dij0,  pawtab(ii)%dij0)
3563  !    end if
3564  !   end do
3565  ! end if
3566 
3567  call ddb_hdr_free(ddb_hdr)
3568 
3569 end subroutine ddb_to_dtset

m_ddb/ddb_type [ Types ]

[ Top ] [ m_ddb ] [ Types ]

NAME

 ddb_type

FUNCTION

  Provides methods to extract and postoprocess the results in the derivative database (DDB)

SOURCE

 83  type,public :: ddb_type
 84 
 85   integer :: msize
 86   ! Maximum size of dynamical matrices and other perturbations (ddk, dde...)
 87 
 88   integer :: mpert
 89   ! TODO: Write function that returns mpert from natom!
 90 
 91   integer :: nblok
 92   ! Number of 2dte blocks in present object
 93 
 94   integer :: natom
 95   ! Number of atoms in the unit cell.
 96 
 97   integer :: ntypat
 98   ! Number of type of atoms.
 99 
100   integer :: occopt
101   ! Occupation option.
102 
103   integer :: prtvol
104   ! Verbosity level.
105 
106   ! These values are used to call the anaddb routines that don't use rprimd, gprimd.
107   real(dp) :: rprim(3,3)
108   real(dp) :: gprim(3,3)
109   real(dp) :: acell(3)
110 
111   integer,allocatable :: flg(:,:)
112   ! flg(msize,nblok)
113   ! flag to indicate presence of a given block
114 
115   integer,allocatable :: typ(:)
116   ! typ(nblok)
117   ! type of each block - ddk, dde, phonon etc...
118 
119   real(dp),allocatable :: amu(:)
120   ! amu(ntypat)
121   ! mass of the atoms (atomic mass unit)
122 
123   real(dp),allocatable :: nrm(:,:)
124   ! nrm(3,nblok)
125   ! norm of the q-points for each block - can be 0 to indicate a direction of approach to gamma
126 
127   real(dp),allocatable :: qpt(:,:)
128   ! qpt(9,nblok)
129   ! q-point vector in reciprocal space (reduced lattice coordinates) for each block
130 
131   real(dp),allocatable :: val(:,:,:)
132   ! val(2,msize,nblok)
133   ! values of the second energy derivatives in each block
134 
135  end type ddb_type
136 
137  public :: ddb_from_file            ! Construct the object from the DDB file.
138  public :: ddb_free                 ! Free dynamic memory.
139  public :: ddb_malloc               ! Allocate dynamic memory
140  public :: ddb_bcast                ! Broadcast the object.
141  public :: ddb_copy                 ! Copy the object.
142  public :: ddb_get_etotal           ! Read the GS total energy.
143  public :: ddb_get_dielt_zeff       ! Reads the Dielectric Tensor and the Effective Charges
144  public :: ddb_get_dielt            ! Reads the Dielectric Tensor
145  public :: ddb_get_dchidet          ! Reads the non-linear optical susceptibility tensor and the
146                                     ! first-order change in the linear dielectric susceptibility
147  public :: ddb_diagoq               ! Compute the phonon frequencies at the specified q-point by performing
148                                     ! a direct diagonalizatin of the dynamical matrix.
149  public :: ddb_get_asrq0            ! Return object used to enforce the acoustic sum rule
150                                     ! from the Dynamical matrix at Gamma. Used in ddb_diagoq.
151  public :: ddb_to_dtset             ! Transfer ddb_hdr to dtset datatype
152 
153  public :: mblktyp1                 ! This routine merges the derivative databases of type 0-4:
154 
155  ! TODO: This routine is deprecated and will be removed
156  public :: mblktyp5                 ! This routine merges the derivative databases of type 5:
157 
158  ! TODO: Add option to change amu.
159  !public :: ddb_change_amu
160  !public :: ddb_print

m_ddb/ddb_write_blok [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 ddb_write_blok

FUNCTION

 This routine writes blocks of data in the DDBs.

INPUTS

 choice= (2 => write), (3 => write minimal info )
 mpert =maximum number of ipert
 msize=maximum size of the arrays flags and values
 nunit=unit number for the data block file

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output
 ddb = ddb block datastructure
 ddb%typ=type of the block:
   0 => total energy
   1 => second-order energy derivatives, non-stationary block
   2 => second-order energy derivatives, stationary block
   3 => third-order energy derivatives
   4 => first-order energy derivatives: forces, stresses and polarization
   5 => second-order eigenvalue derivatives
 ddb%flg(msize)=flag for every matrix element (0=> the element is
  not in the data block), (1=> the element is in the data blok)
 ddb%qpt(9)=wavevector of the perturbation(s). The elements from
  1 to 3 are used if we are dealing with the 2nd derivative of
  total energy (only one wavevector), while all elements are
  used in case of a third order derivative of total energy
  (three wavevector could be present)
 ddb%nrm(3)=normalization factors for the three allowed wavevectors.
 ddb%val(2,msize)=real(dp), complex, value of the
  matrix elements that are present in the data block
 blkval2(2,msize,mband,nkpt) = value of the matrix elements
  that are present in a block of EIGR2D/EIGI2D

NOTES

 only executed by one processor.

PARENTS

      ddb_interpolate,dfptnl_doutput,gstate,mblktyp1,mblktyp5

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

3104 subroutine ddb_write_blok(ddb,iblok,choice,mband,mpert,msize,nkpt,nunit,&
3105 &     blkval2,kpt) !optional
3106 
3107 
3108 !This section has been created automatically by the script Abilint (TD).
3109 !Do not modify the following lines by hand.
3110 #undef ABI_FUNC
3111 #define ABI_FUNC 'ddb_write_blok'
3112 !End of the abilint section
3113 
3114  implicit none
3115 
3116 !Arguments -------------------------------
3117 !scalars
3118  integer,intent(in) :: choice,mband,mpert,msize,nkpt,nunit
3119  integer,intent(in) :: iblok
3120  type(ddb_type),intent(in) :: ddb
3121 !arrays
3122  real(dp),intent(in),optional :: kpt(3,nkpt)
3123  real(dp),intent(in),optional :: blkval2(2,msize,mband,nkpt)
3124 
3125 !Local variables -------------------------
3126 !scalars
3127  integer :: iband,idir1,idir2,idir3,ii,ikpt,ipert1,ipert2,ipert3
3128  integer :: nelmts
3129 
3130 ! *********************************************************************
3131 
3132 
3133 !Count the number of elements
3134  nelmts=0
3135  do ii=1,msize
3136    if(ddb%flg(ii,iblok)==1)nelmts=nelmts+1
3137  end do
3138 
3139 !Write the block type and number of elements
3140  write(nunit,*)' '
3141  if (ddb%typ(iblok) == 0) then
3142    write(nunit, '(a,i8)' )' Total energy                 - # elements :',nelmts
3143  else if (ddb%typ(iblok)==1) then
3144    write(nunit, '(a,i8)' )' 2nd derivatives (non-stat.)  - # elements :',nelmts
3145  else if(ddb%typ(iblok)==2) then
3146    write(nunit, '(a,i8)' )' 2nd derivatives (stationary) - # elements :',nelmts
3147  else if(ddb%typ(iblok)==3) then
3148    write(nunit, '(a,i8)' )' 3rd derivatives              - # elements :',nelmts
3149  else if (ddb%typ(iblok) == 4) then
3150    write(nunit, '(a,i8)' )' 1st derivatives              - # elements :',nelmts
3151  else if (ddb%typ(iblok) == 5) then
3152    write(nunit, '(a,i8)' )' 2nd eigenvalue derivatives   - # elements :',nelmts
3153  end if
3154 
3155 !Write the 2nd derivative block
3156  if(ddb%typ(iblok)==1.or.ddb%typ(iblok)==2)then
3157 
3158 !  Write the phonon wavevector
3159    write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
3160 
3161 !  Write the matrix elements
3162    if(choice==2)then
3163      ii=0
3164      do ipert2=1,mpert
3165        do idir2=1,3
3166          do ipert1=1,mpert
3167            do idir1=1,3
3168              ii=ii+1
3169              if(ddb%flg(ii,iblok)==1)then
3170                write(nunit,'(4i4,2d22.14)')idir1,ipert1,idir2,ipert2,ddb%val(1,ii,iblok),ddb%val(2,ii,iblok)
3171              end if
3172            end do
3173          end do
3174        end do
3175      end do
3176    end if
3177 
3178 !  Write the 3rd derivative block
3179  else if(ddb%typ(iblok)==3)then
3180 
3181 !  Write the phonon wavevectors
3182    write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
3183    write(nunit, '(a,3es16.8,f6.1)' )'    ',(ddb%qpt(ii,iblok),ii=4,6),ddb%nrm(2,iblok)
3184    write(nunit, '(a,3es16.8,f6.1)' )'    ',(ddb%qpt(ii,iblok),ii=7,9),ddb%nrm(3,iblok)
3185 
3186 !  Write the matrix elements
3187    if(choice==2)then
3188      ii=0
3189      do ipert3=1,mpert
3190        do idir3=1,3
3191          do ipert2=1,mpert
3192            do idir2=1,3
3193              do ipert1=1,mpert
3194                do idir1=1,3
3195                  ii=ii+1
3196                  if(ddb%flg(ii,iblok)==1)then
3197                    write(nunit, '(6i4,2d22.14)' )&
3198 &                   idir1,ipert1,idir2,ipert2,idir3,ipert3,ddb%val(1,ii,iblok),ddb%val(2,ii,iblok)
3199                  end if
3200                end do
3201              end do
3202            end do
3203          end do
3204        end do
3205      end do
3206    end if
3207 
3208 !  Write total energy
3209  else if (ddb%typ(iblok) == 0) then
3210    if (choice == 2) then
3211      write(nunit,'(2d22.14)')ddb%val(1,1,iblok),ddb%val(2,1,iblok)
3212    end if
3213 
3214 !  Write the 1st derivative blok
3215  else if (ddb%typ(iblok) == 4) then
3216    if (choice == 2) then
3217      ii = 0
3218      do ipert1 = 1, mpert
3219        do idir1 = 1, 3
3220          ii = ii + 1
3221          if (ddb%flg(ii,iblok) == 1) then
3222            write(nunit,'(2i4,2d22.14)')idir1,ipert1,ddb%val(1,ii,iblok),ddb%val(2,ii,iblok)
3223          end if
3224        end do
3225      end do
3226    end if
3227 
3228  else if (ddb%typ(iblok)==5) then
3229 !  Write the phonon wavevector
3230    write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok)
3231 !  Write the matrix elements
3232    if(choice==2)then
3233      if(present(blkval2).and.present(kpt))then
3234        do ikpt=1,nkpt
3235          write(nunit,'(a,3es16.8)')' K-point:',(kpt(ii,ikpt),ii=1,3)
3236          do iband=1,mband
3237            write(nunit,'(a,i3)')' Band:',iband
3238            ii=0
3239            do ipert2=1,mpert
3240              do idir2=1,3
3241                do ipert1=1,mpert
3242                  do idir1=1,3
3243                    ii=ii+1
3244                    if(ddb%flg(ii,iblok)==1)then
3245                      write(nunit,'(4i4,2d22.14)')idir1,ipert1,idir2,ipert2,&
3246 &                     blkval2(1,ii,iband,ikpt),blkval2(2,ii,iband,ikpt)
3247                    end if
3248                  end do !idir1
3249                end do  !ipert1
3250              end do   !idir2
3251            end do    !ipert2
3252          end do     !iband
3253        end do      !ikpt
3254      end if !blkval2
3255    end if !choice
3256  end if !ddb%typ(iblok)
3257 
3258 end subroutine ddb_write_blok

m_ddb/dfptnl_doutput [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 dfptnl_doutput

FUNCTION

 Write the matrix of third-order derivatives to the output file and the DDB

INPUTS

  blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte
   has been calculated ; 0 otherwise )
  d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
  mpert =maximum number of ipert
  natom=Number of atoms
  ntypat=Number of type of atoms
  unddb = unit number for DDB output

NOTES

  d3 holds the third-order derivatives before computing
  the permutations of the perturbations.

PARENTS

      nonlinear

CHILDREN

      ddb_free,ddb_malloc,ddb_write_blok,wrtout

SOURCE

3289 subroutine dfptnl_doutput(blkflg,d3,mband,mpert,nkpt,natom,ntypat,unddb)
3290 
3291 
3292 !This section has been created automatically by the script Abilint (TD).
3293 !Do not modify the following lines by hand.
3294 #undef ABI_FUNC
3295 #define ABI_FUNC 'dfptnl_doutput'
3296 !End of the abilint section
3297 
3298  implicit none
3299 
3300 !Arguments -------------------------------
3301 !scalars
3302  integer,intent(in) :: mband,mpert,nkpt,unddb,natom,ntypat
3303 !arrays
3304  integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert)
3305  real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert)
3306 
3307 !Local variables -------------------------
3308 !scalars
3309  integer :: choice,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,index,msize
3310  character(len=500) :: message
3311  type(ddb_type) :: ddb
3312 
3313 !*************************************************************************
3314 
3315  msize = 27*mpert*mpert*mpert
3316  call ddb_malloc(ddb,msize,1,natom,ntypat)
3317 
3318  choice = 2
3319 
3320  ddb%typ = 3
3321  ddb%nrm = one
3322  ddb%qpt = zero   ! this has to be changed in case anharmonic
3323 !force constants have been computed
3324 
3325 
3326 !Write blok of third-order derivatives to ouput file
3327 
3328  write(message,'(a,a,a,a,a)')ch10,&
3329 & ' Matrix of third-order derivatives (reduced coordinates)',ch10,&
3330 & ' before computing the permutations of the perturbations',ch10
3331  call wrtout(ab_out,message,'COLL')
3332 
3333  write(ab_out,*)'    j1       j2       j3              matrix element'
3334  write(ab_out,*)' dir pert dir pert dir pert           real part           imaginary part'
3335 
3336  do i1pert=1,mpert
3337    do i1dir=1,3
3338      do i2pert=1,mpert
3339        do i2dir=1,3
3340          do i3pert=1,mpert
3341            do i3dir=1,3
3342 
3343              index = i1dir + &
3344 &             3*((i1pert-1)+mpert*((i2dir-1) + &
3345 &             3*((i2pert-1)+mpert*((i3dir-1) + 3*(i3pert-1)))))
3346              ddb%flg(index,1) = blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3347              ddb%val(:,index,1)= d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3348 
3349              if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=0) then
3350 
3351                write(ab_out,'(3(i4,i5),2f22.10)')&
3352 &               i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,&
3353 &               d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3354 
3355              end if
3356 
3357            end do
3358          end do
3359        end do
3360      end do
3361    end do
3362  end do
3363 
3364 !Write blok of third-order derivatives to DDB
3365  call ddb_write_blok(ddb,1,choice,mband,mpert,msize,nkpt,unddb)
3366 
3367  call ddb_free(ddb)
3368 
3369 end subroutine dfptnl_doutput

m_ddb/dtchi [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 dtchi

FUNCTION

 Reads the non-linear optical susceptibility tensor and the
 first-order change in the linear dielectric susceptibility
 induced by an atomic displacement in the Gamma Block coming from the Derivative Data Base
 (third-order derivatives).

INPUTS

 blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies
 natom= number of atoms in unit cell
 mpert =maximum number of ipert
 ramansr= if /= 0, impose sum rule on first-order derivatives
                   of the electronic susceptibility with respect
                   to atomic displacements
 nlflag= if =3, only the non-linear optical susceptibilities is computed

OUTPUT

 dchide(3,3,3) = non-linear optical coefficients
 dchidt(natom,3,3,3) = first-order change of the electronic dielectric
   tensor induced by an individual atomic displacement

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

2190 subroutine dtchi(blkval,dchide,dchidt,mpert,natom,ramansr,nlflag)
2191 
2192 
2193 !This section has been created automatically by the script Abilint (TD).
2194 !Do not modify the following lines by hand.
2195 #undef ABI_FUNC
2196 #define ABI_FUNC 'dtchi'
2197 !End of the abilint section
2198 
2199  implicit none
2200 
2201 !Arguments -------------------------------
2202 !scalars
2203  integer,intent(in) :: mpert,natom,ramansr,nlflag
2204 !arrays
2205  real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert)
2206  real(dp),intent(out) :: dchide(3,3,3),dchidt(natom,3,3,3)
2207 
2208 !Local variables -------------------------
2209 !scalars
2210  integer :: depl,elfd1,elfd2,elfd3,iatom,ivoigt
2211  logical :: iwrite
2212  real(dp) :: wttot
2213 !arrays
2214  integer :: voigtindex(6,2)
2215  real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert),dvoigt(3,6),sumrule(3,3,3)
2216  real(dp) :: wghtat(natom)
2217 
2218 ! *********************************************************************
2219 
2220  d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/))
2221  d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/))
2222 
2223 !Extraction of non-linear optical coefficients
2224  do elfd1 = 1,3
2225    do elfd2 = 1,3
2226      do elfd3 = 1,3
2227        dchide(elfd1,elfd2,elfd3) = d3cart(1,elfd1,natom+2,elfd2,natom+2,elfd3,natom+2)
2228      end do
2229    end do
2230  end do
2231 
2232 !Transform to Voigt notations
2233  voigtindex(:,1) = (/1,2,3,2,1,1/)
2234  voigtindex(:,2) = (/1,2,3,3,3,2/)
2235  do ivoigt = 1, 6
2236    elfd2 = voigtindex(ivoigt,1)
2237    elfd3 = voigtindex(ivoigt,2)
2238    do elfd1 = 1, 3
2239      dvoigt(elfd1,ivoigt) = 0.5_dp*(dchide(elfd1,elfd2,elfd3) + dchide(elfd1,elfd3,elfd2))
2240    end do
2241  end do
2242 
2243 !Transform to pm/V
2244  dvoigt(:,:) = dvoigt(:,:)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
2245 
2246 !Extraction of $\frac{d \chi}{d \tau}$
2247  if (nlflag < 3) then
2248    do iatom = 1, natom
2249      do depl = 1,3
2250        do elfd1 = 1,3
2251          do elfd2 = 1,3
2252            dchidt(iatom,depl,elfd1,elfd2) = d3cart(1,depl,iatom,elfd1,natom+2,elfd2,natom+2)
2253          end do
2254        end do
2255      end do
2256    end do
2257  end if
2258 
2259  wghtat(:) = zero
2260  if (ramansr == 1) then
2261    wghtat(:) = one/dble(natom)
2262 
2263  else if (ramansr == 2) then
2264 
2265    wttot = zero
2266    do iatom = 1, natom
2267      do depl = 1,3
2268        do elfd1 = 1,3
2269          do elfd2 = 1,3
2270            wghtat(iatom) = wghtat(iatom) + abs(dchidt(iatom,depl,elfd1,elfd2))
2271          end do
2272        end do
2273      end do
2274      wttot = wttot + wghtat(iatom)
2275    end do
2276 
2277    wghtat(:) = wghtat(:)/wttot
2278  end if
2279 
2280  iwrite = ab_out > 0
2281 
2282  if (iwrite) then
2283    write(ab_out,*)ch10
2284    write(ab_out,*)'Non-linear optical coefficients d (pm/V)'
2285    write(ab_out,'(6f12.6)')dvoigt(1,:)
2286    write(ab_out,'(6f12.6)')dvoigt(2,:)
2287    write(ab_out,'(6f12.6)')dvoigt(3,:)
2288  end if
2289 
2290  if (ramansr /= 0) then
2291    if (iwrite) then
2292      write(ab_out,*)ch10
2293      write(ab_out,*)'The violation of the Raman sum rule'
2294      write(ab_out,*)'by the first-order electronic dielectric tensors ','is as follows'
2295      write(ab_out,*)'    atom'
2296      write(ab_out,*)' displacement'
2297    end if
2298 
2299    sumrule(:,:,:) = zero
2300    do elfd2 = 1,3
2301      do elfd1 = 1,3
2302        do depl = 1,3
2303          do iatom = 1, natom
2304            sumrule(depl,elfd1,elfd2) = sumrule(depl,elfd1,elfd2) + &
2305 &           dchidt(iatom,depl,elfd1,elfd2)
2306          end do
2307          do iatom = 1, natom
2308            dchidt(iatom,depl,elfd1,elfd2) = dchidt(iatom,depl,elfd1,elfd2) - &
2309 &           wghtat(iatom)*sumrule(depl,elfd1,elfd2)
2310          end do
2311        end do
2312      end do
2313    end do
2314 
2315    if (iwrite) then
2316      do depl = 1,3
2317        write(ab_out,'(6x,i2,3(3x,f16.9))') depl,sumrule(depl,1,1:3)
2318        write(ab_out,'(8x,3(3x,f16.9))') sumrule(depl,2,1:3)
2319        write(ab_out,'(8x,3(3x,f16.9))') sumrule(depl,3,1:3)
2320        write(ab_out,*)
2321      end do
2322     end if
2323  end if    ! ramansr
2324 
2325  if (nlflag < 3) then
2326    if (iwrite) then
2327      write(ab_out,*)ch10
2328      write(ab_out,*)' First-order change in the electronic dielectric '
2329      write(ab_out,*)' susceptibility tensor (Bohr^-1)'
2330      write(ab_out,*)' induced by an atomic displacement'
2331      if (ramansr /= 0) then
2332        write(ab_out,*)' (after imposing the sum over all atoms to vanish)'
2333      end if
2334      write(ab_out,*)'  atom  displacement'
2335 
2336      do iatom = 1,natom
2337        do depl = 1,3
2338          write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')iatom,depl,dchidt(iatom,depl,1,:)
2339          write(ab_out,'(16x,3(3x,f16.9))')dchidt(iatom,depl,2,:)
2340          write(ab_out,'(16x,3(3x,f16.9))')dchidt(iatom,depl,3,:)
2341        end do
2342 
2343        write(ab_out,*)
2344      end do
2345    end if
2346  end if
2347 
2348 end subroutine dtchi

m_ddb/dtech9 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 dtech9

FUNCTION

 Reads the Dielectric Tensor and the Effective Charges in the
 Gamma Block coming from the Derivative Data Base.

INPUTS

 natom= number of atoms in unit cell
 iblok= index of the Gamma block
 mpert =maximum number of ipert
 nblok= number of blocks in the DDB
 blkval(2,3*mpert*3*mpert,nblok)=  dynamical matrices
  In our case, the nblok is restricted to iblok

OUTPUT

 zeff(3,3,natom)=effective charge on each atom, versus electric
  field and atomic displacement. Note the following convention:
  zeff(electric field direction, atomic direction, atom number)
 dielt(3,3)=dielectric tensor

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

2091 subroutine dtech9(blkval,dielt,iblok,mpert,natom,nblok,zeff)
2092 
2093 
2094 !This section has been created automatically by the script Abilint (TD).
2095 !Do not modify the following lines by hand.
2096 #undef ABI_FUNC
2097 #define ABI_FUNC 'dtech9'
2098 !End of the abilint section
2099 
2100  implicit none
2101 
2102 !Arguments -------------------------------
2103 !scalars
2104  integer,intent(in) :: iblok,mpert,natom,nblok
2105 !arrays
2106  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok)
2107  real(dp),intent(out) :: dielt(3,3),zeff(3,3,natom)
2108 
2109 !Local variables -------------------------
2110 !scalars
2111  integer :: depl,elec,elec1,elec2,iatom
2112  character(len=1000) :: message
2113 
2114 ! *********************************************************************
2115 
2116 !Extration of effectives charges
2117  do iatom=1,natom
2118    do elec=1,3
2119      do depl=1,3
2120        zeff(elec,depl,iatom)=0.5*&
2121 &       (blkval(1,depl,iatom,elec,natom+2,iblok)+&
2122 &       blkval(1,elec,natom+2,depl,iatom,iblok))
2123      end do
2124    end do
2125  end do
2126 
2127 !Extration of dielectric tensor
2128  do elec1=1,3
2129    do elec2=1,3
2130      dielt(elec1,elec2)=blkval(1,elec1,natom+2,elec2,natom+2,iblok)
2131    end do
2132  end do
2133 
2134  write(message,'(a,3es16.6,3es16.6,3es16.6)' )&
2135 & ' Dielectric Tensor ',&
2136 & dielt(1,1),dielt(1,2),dielt(1,3),&
2137 & dielt(2,1),dielt(2,2),dielt(2,3),&
2138 & dielt(3,1),dielt(3,2),dielt(3,3)
2139 
2140  call wrtout(std_out,message,'COLL')
2141 
2142  write(message,'(a)' ) ' Effectives Charges '
2143  call wrtout(std_out,message,'COLL')
2144  do iatom=1,natom
2145    write(message,'(a,i4,3es16.6,3es16.6,3es16.6)' )' atom ',iatom,&
2146 &   zeff(1,1,iatom),zeff(1,2,iatom),zeff(1,3,iatom),&
2147 &   zeff(2,1,iatom),zeff(2,2,iatom),zeff(2,3,iatom),&
2148 &   zeff(3,1,iatom),zeff(3,2,iatom),zeff(3,3,iatom)
2149     call wrtout(std_out,message,'COLL')
2150  end do
2151 
2152 end subroutine dtech9

m_ddb/gamma9 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 gamma9

FUNCTION

 This small routine checks if the wavevector qphon and the
 corresponding normalisation factor represent a phonon at Gamma.

INPUTS

 qphon(3)=wavevector
 qphnrm=normalisation factor
 qtol=tolerance

OUTPUT

 gamma= if 1, means that the wavevector is indeed at Gamma otherwise 0.

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

799 subroutine gamma9(gamma,qphon,qphnrm,qtol)
800 
801 
802 !This section has been created automatically by the script Abilint (TD).
803 !Do not modify the following lines by hand.
804 #undef ABI_FUNC
805 #define ABI_FUNC 'gamma9'
806 !End of the abilint section
807 
808  implicit none
809 
810 !Arguments -------------------------------
811 !scalars
812  integer,intent(out) :: gamma
813  real(dp),intent(in) :: qphnrm,qtol
814 !arrays
815  real(dp),intent(in) :: qphon(3)
816 
817 ! *********************************************************************
818 
819  if( (abs(qphon(1))<qtol .and. abs(qphon(2))<qtol .and. abs(qphon(3))<qtol) .or. abs(qphnrm)<qtol ) then
820    gamma=1
821  else
822    gamma=0
823  end if
824 
825 end subroutine gamma9

m_ddb/gtblk9 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 gtblk9

FUNCTION

 This routine (get block) finds the block that contains the
 information on the derivatives of the total energy specified
 by the parameters rfphon,rfelfd,rfstrs,rftyp and
 the phonon wavevectors qphon (and their normalisation).
 In case the DDB does not contain this information, the subroutine returns iblok=0

INPUTS

 ddb = ddb blok datastructure
   flg(msize,nblok)=flag for every matrix element:
      0 => the element is not in the data block.
      1 => the element is in the data blok.
   nrm(3,nblok)=normalization factors for the three allowed wavevectors
   qpt(3,nblok)=wavevector of the perturbation(s). The elements
   typ(nblok)=type of the block.
      (1=> non-stationary block),
      (2=> stationary block),
      (3=> third order derivative).
 qphon(3,3)=wavevectors for the three possible phonons
  (note : only one should be used in case of second derivative of total energy,
  because we know that the second is the opposite of this value)
 qphnrm(3) =normalisation factors for the three possible phonons
 rfphon(4) = 1=> response to phonons (for the four possible derivatives. Two should be used for a second derivative of total energy)
 rfelfd(4) = 1=> d/dk, 2=> electric field only, 3=> both (see comment on rfphon)
 rfstrs(4) = 1=> uniaxial stresses, 2=> shear stresses, 3=> both (see comment on rfphon)
 rftyp =
   0 => total energy
   1 => non-stationary formulation of the 2nd derivative
   2 => stationary formulation of the 2nd derivative
   3 => third derivative of total energy
   4 => first-order derivatives of total energy

OUTPUT

 iblok= number of the block that corresponds to the specifications

PARENTS

      anaddb,ddb_interpolate,m_ddb,m_effective_potential_file,m_phonons
      thmeig

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

533 subroutine gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
534 
535 
536 !This section has been created automatically by the script Abilint (TD).
537 !Do not modify the following lines by hand.
538 #undef ABI_FUNC
539 #define ABI_FUNC 'gtblk9'
540 !End of the abilint section
541 
542  implicit none
543 
544 !Arguments -------------------------------
545 !scalars
546  integer,intent(in) :: rftyp
547  integer,intent(out) :: iblok
548  type(ddb_type),intent(in) :: ddb
549 !arrays
550  integer,intent(in) :: rfelfd(4),rfphon(4),rfstrs(4)
551  real(dp),intent(inout) :: qphnrm(3),qphon(3,3)
552 
553 !Local variables -------------------------
554 !scalars
555  integer :: blkgam,ider,idir,idir1,idir2,idir3,ii,index,ipert,ipert1,ipert2
556  integer :: ipert3,nder,ok,mpert,natom
557  character(len=500) :: message
558 !arrays
559  integer :: gamma(3)
560  integer,allocatable :: worki(:,:)
561  real(dp) :: qpt(3)
562 
563 ! *********************************************************************
564 
565  mpert = ddb%mpert
566  natom = ddb%natom
567 
568 !Get the number of derivative
569  if(rftyp==1.or.rftyp==2)then
570    nder=2
571  else if(rftyp==3)then
572    nder=3
573  else if(rftyp==0)then
574    nder=0
575  else if(rftyp==4)then
576    nder=1
577  else
578    write(message, '(a,i0,a)')' rftyp is equal to ',rftyp,'. The only allowed values are 0, 1, 2, 3 or 4.'
579    MSG_BUG(message)
580  end if
581 
582 !In case of a second-derivative, a second phonon wavevector is provided.
583  if(nder==2)then
584    do ii=1,3
585      qphon(ii,2)=-qphon(ii,1)
586    end do
587    qphnrm(2)=qphnrm(1)
588  end if
589 
590 !In case of a third derivative, the sum of wavevectors to gamma is checked
591  if (nder == 3) then
592    qpt(:) = qphon(:,1)/qphnrm(1) + qphon(:,2)/qphnrm(2) + qphon(:,3)/qphnrm(3)
593    call gamma9(gamma(nder),qpt,qphnrm(1),DDB_QTOL)
594    if (gamma(nder) == 0) then
595      write(message,'(a,a,a)')&
596 &     'the sum of the wavevectors of the third-order energy is ',ch10,&
597 &     'not equal to zero'
598      MSG_ERROR(message)
599    end if
600  end if
601 
602 !Check the validity of the requirement
603  do ider=1,nder
604    ! Identifies if qphon is at gamma
605    call gamma9(gamma(ider),qphon(1:3,ider),qphnrm(ider),DDB_QTOL)
606 
607    if(gamma(ider)==0)then
608      if(rfstrs(ider)/=0.or.rfelfd(ider)/=0)then
609        write(message, '(a,a)' )&
610 &       'Not yet able to handle stresses or electric fields',ch10,&
611 &       'with non-zero wavevector.'
612        MSG_BUG(message)
613      end if
614    end if
615  end do
616 
617 !Initialise the perturbation table
618  ABI_MALLOC(worki,(mpert,4))
619  worki(:,1:nder)=0
620 
621 !Build the perturbation table
622  do ider=1,nder
623 !  First the phonons
624    if(rfphon(ider)==1)then
625      do ipert=1,natom
626        worki(ipert,ider)=1
627      end do
628    end if
629 !  Then the d/dk
630    if(rfelfd(ider)==1.or.rfelfd(ider)==3)then
631      worki(natom+1,ider)=1
632    end if
633 !  Then the electric field
634    if(rfelfd(ider)==2.or.rfelfd(ider)==3)then
635      worki(natom+2,ider)=1
636    end if
637 !  Then the uniaxial stress
638    if(rfstrs(ider)==1.or.rfstrs(ider)==3)then
639      worki(natom+3,ider)=1
640    end if
641 !  At last, the shear stress
642    if(rfstrs(ider)==2.or.rfstrs(ider)==3)then
643      worki(natom+4,ider)=1
644    end if
645  end do
646 
647 !Examine every blok :
648  do iblok=1,ddb%nblok
649 
650 !  If this variable is still 1 at the end of the examination, the blok is the good one...
651    ok=1
652 
653 !  Check the type
654    if(rftyp/=ddb%typ(iblok)) ok=0
655 
656 !  Check the wavevector
657    if( ok==1 )then
658 
659      if (nder == 2) then
660        call gamma9(blkgam,ddb%qpt(1:3,iblok),ddb%nrm(1,iblok),DDB_QTOL)
661        if(blkgam/=gamma(1))then
662          ok=0
663        else if(blkgam==0)then
664          do idir=1,3
665            if( abs( ddb%qpt(idir,iblok)/ddb%nrm(1,iblok) - qphon(idir,1)/qphnrm(1) )>DDB_QTOL ) ok=0
666          end do
667        end if
668 
669      else if (nder == 3) then
670        do ider = 1, nder
671          do idir=1,3
672            if( abs( ddb%qpt(idir+3*(ider-1),iblok)/ddb%nrm(ider,iblok) - &
673 &           qphon(idir,ider)/qphnrm(ider) )>DDB_QTOL )then
674              ok=0
675            end if ! qphon
676          end do ! idir
677        end do ! nder
678      end if  ! nder
679 
680    end if ! ok
681 
682 !  Check if there is enough information in this blok
683    if( ok==1 )then
684 
685      if (nder == 0) then
686        if (ddb%flg(1,iblok) /= 1) then
687          ok = 0
688          if (ddb%prtvol > 1) then
689            write(message,'(a,i0,3a)' )&
690            'The block ',iblok,' does not match the requirement',ch10,&
691            'because it lacks the total energy'
692            MSG_COMMENT(message)
693          end if
694        end if
695      end if
696 
697      do ipert1=1,mpert
698 
699        if ((nder == 4).and.(worki(ipert1,4) == 1).and.(ok == 1)) then
700          do idir1 = 1, 3
701            index = 3*(ipert1 - 1) + idir1
702            if (ddb%flg(index,iblok) /= 1) ok = 0
703          end do
704        end if
705 
706        if (worki(ipert1,1)==1 .and. ok==1 )then
707          do ipert2=1,mpert
708            if (worki(ipert2,2)==1 .and. ok==1 )then
709              do idir1=1,3
710                do idir2=1,3
711 
712                  if (nder == 2) then
713                    index=idir1+ 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
714                    if (ddb%flg(index,iblok)/=1) ok=0
715 
716                  else if (nder == 3) then
717                    do ipert3 = 1, mpert
718                      if (worki(ipert3,3) == 1 .and. ok == 1) then
719                        do idir3 = 1, 3
720                          index = idir1 + &
721 &                         3*((ipert1 - 1) + mpert*((idir2 - 1) + &
722 &                         3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
723                          if (ddb%flg(index,iblok) /= 1) ok = 0
724                        end do  ! idir3
725                      end if ! worki(ipert3,3)
726                    end do ! i3pert
727                  end if
728 
729                end do
730              end do
731            end if
732          end do
733        end if
734      end do
735    end if
736 
737 !  Now that everything has been checked, eventually end the search
738    if(ok==1)exit
739  end do
740 
741  if(ok==0)then
742    iblok=0
743 
744    if (ddb%prtvol > 1) then
745      write(message, '(3a)' )&
746 &     ' gtblk9 : ',ch10,&
747 &     '  Unable to find block corresponding to the following specifications :'
748      call wrtout(std_out,message,'COLL')
749      write(message, '(a,i3)' )' Type (rfmeth) =',rftyp
750      call wrtout(std_out,message,'COLL')
751      write(message, '(a)' ) ' ider qphon(3)         qphnrm   rfphon rfelfd rfstrs'
752      call wrtout(std_out,message,'COLL')
753      do ider=1,nder
754        write(message, '(i4,4f6.2,3i7)' )&
755        ider,(qphon(ii,ider),ii=1,3),qphnrm(ider),rfphon(ider),rfelfd(ider),rfstrs(ider)
756        call wrtout(std_out,message,'COLL')
757      end do
758    end if
759  end if
760 
761  if (ok==1 .and. ddb%prtvol > 1) then
762    write(message,'(a,i0,2a)')' gtblk9: found block number ',iblok,' agree with',' specifications '
763    call wrtout(std_out,message,'COLL')
764  end if
765 
766  ABI_FREE(worki)
767 
768 end subroutine gtblk9

m_ddb/mblktyp1 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 mblktyp1

FUNCTION

 This routine merges the derivative databases of type 0-4:
 Total energy, (2nd derivatives (non-stat.),2nd derivatives (stationary),
 3rd derivatives, 1st derivatives

 The heading of the database is read, then the heading
 of the temporary database to be added is read,
 the code check their compatibility, and create a new
 database that mixes the old and the temporary ones.
 This process can be iterated.
 The whole database will be stored in central memory.

INPUTS

     chkopt=option for consistency checks between DDB files
     ddbun=define input and output unit numbers
     dscrpt=description of the output file
     filnam=name of input or output file
     mddb=maximum number of databases (cannot be made dynamic)
     nddb=number of input DDBs
     vrsddb=current version of the DDB

OUTPUT

     msym=maximum number of symmetry elements in space group
     Merge the file

PARENTS

      mrgddb

CHILDREN

      ddb_free,ddb_hdr_compare,ddb_hdr_free,ddb_hdr_open_read
      ddb_hdr_open_write,ddb_malloc,ddb_write_blok,read_blok8,wrtout

SOURCE

3613 subroutine mblktyp1(chkopt,ddbun,dscrpt,filnam,mddb,msym,nddb,vrsddb)
3614 
3615 
3616 !This section has been created automatically by the script Abilint (TD).
3617 !Do not modify the following lines by hand.
3618 #undef ABI_FUNC
3619 #define ABI_FUNC 'mblktyp1'
3620 !End of the abilint section
3621 
3622  implicit none
3623 
3624 !Arguments -------------------------------
3625 !scalars
3626  integer,intent(in) :: chkopt,ddbun,mddb,nddb,vrsddb
3627  integer,intent(out) :: msym
3628  character(len=fnlen),intent(in) :: dscrpt
3629  character(len=fnlen),intent(in) :: filnam(mddb+1)
3630 
3631 !Local variables -------------------------
3632 !scalars
3633  integer :: choice,dimekb,iblok,iblok1,iblok2
3634  integer :: iddb,ii,lmnmax,matom
3635  integer :: mband,mblktyp,mblok,mkpt,mpert,msize,mtypat
3636  integer :: nblok,nblokt,nq
3637  integer :: tmerge,usepaw
3638  real(dp),parameter :: qtol=2.0d-8
3639  real(dp) :: diff
3640  type(ddb_type) :: ddb
3641  type(ddb_hdr_type) :: ddb_hdr, ddb_hdr8
3642  character(len=500) :: message
3643 !arrays
3644  integer,allocatable :: mgblok(:)!,lloc(:)
3645 
3646 ! *********************************************************************
3647 
3648  ! Make sure there is more than one ddb to be read
3649  if(nddb==1)then
3650    write(message, '(a,a,a,a,a)' )&
3651 &   'The initialisation mode of MRGDDB, that uses nddb=1,',&
3652 &   'has been disabled in version 2.1 of ABINIT.',&
3653 &   'Action: you should use DDBs that include the symmetry',&
3654 &   'information (and that can be used and merged without',&
3655 &   'initialisation), or you should use ABINITv2.0.'
3656    MSG_ERROR(message)
3657  end if
3658 
3659 !Evaluate the maximal dimensions of arrays
3660  dimekb=0 ; matom=0 ; mband=0  ; mblok=0 ; mkpt=0
3661  msize=0  ; mtypat=0 ; lmnmax=0 ; usepaw=0 ; mblktyp = 1
3662  msym=192
3663 
3664  do iddb=1,nddb
3665    call ddb_hdr_open_read(ddb_hdr, filnam(iddb+1), ddbun, vrsddb, dimonly=1)
3666 
3667    mblok=mblok+ddb_hdr%nblok
3668    mblktyp=max(mblktyp,ddb_hdr%mblktyp)
3669    matom=max(matom,ddb_hdr%matom)
3670    mkpt=max(mkpt,ddb_hdr%mkpt)
3671    mtypat=max(mtypat,ddb_hdr%mtypat)
3672    msym=max(msym,ddb_hdr%msym)
3673    mband=max(mband,ddb_hdr%mband)
3674    dimekb=max(dimekb,ddb_hdr%psps%dimekb)
3675    lmnmax=max(lmnmax,ddb_hdr%psps%lmnmax)
3676    usepaw=max(usepaw,ddb_hdr%usepaw)
3677 
3678    call ddb_hdr_free(ddb_hdr)
3679  end do
3680 
3681  mpert=matom+6
3682  msize=3*mpert*3*mpert
3683  if(mblktyp==3)msize=msize*3*mpert
3684 
3685  call ddb_malloc(ddb,msize,mblok,matom,mtypat)
3686 
3687 !Allocate arrays
3688  ABI_ALLOCATE(mgblok,(mblok))
3689 
3690 !**********************************************************************
3691 
3692 !Read the first database
3693 
3694  write(std_out,*)' read the input derivative database information'
3695  call ddb_hdr_open_read(ddb_hdr, filnam(2), ddbun, vrsddb, &
3696 & matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,&
3697 & msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw)
3698 
3699  if(ddb_hdr%nblok>=1)then
3700 !  Read the blocks from the input database.
3701    write(message, '(a,i5,a)' ) ' read ',ddb_hdr%nblok, ' blocks from the input DDB '
3702    call wrtout(std_out,message,'COLL')
3703    do iblok=1,ddb_hdr%nblok
3704      call read_blok8(ddb,iblok,ddb_hdr%nband(1),mpert,msize,ddb_hdr%nkpt,ddbun)
3705 !    Setup merged indicator
3706      mgblok(iblok)=0
3707    end do
3708  else
3709    write(message, '(a)' )' No bloks in the first ddb '
3710    call wrtout(std_out,message,'COLL')
3711  end if
3712 
3713 !Close the first ddb
3714  close(ddbun)
3715 
3716 !*********************************************
3717 
3718  nblok = ddb_hdr%nblok
3719 !In case of merging of DDBs, iterate the reading
3720  do iddb=2,nddb
3721 
3722 !  Open the corresponding input DDB, and read the database file information
3723    write(message, '(a,a,i6)' )ch10,' read the input derivative database number',iddb
3724    call wrtout(std_out,message,'COLL')
3725 
3726    call ddb_hdr_open_read(ddb_hdr8, filnam(iddb+1), ddbun, vrsddb, &
3727 &   matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,&
3728 &   msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw)
3729 
3730    if (chkopt==1)then
3731 !    Compare the current DDB and input DDB information.
3732 !    In case of an inconsistency, halt the execution.
3733      call wrtout(std_out, ' compare the current and input DDB information', 'COLL')
3734      call ddb_hdr_compare(ddb_hdr, ddb_hdr8)
3735 
3736    else if(chkopt==0)then
3737 !    No comparison between the current DDB and input DDB information.
3738      write(message, '(a)' )' no comparison between the current and input DDB information'
3739      call wrtout(std_out,message,'COLL')
3740      write(message, '(5a)' )&
3741 &     'No comparison/check is performed for the current and input DDB information ',ch10,&
3742 &     'because argument --nostrict was passed to the command line. ',ch10,&
3743 &     'Use at your own risk!'
3744      MSG_COMMENT(message)
3745    end if
3746 
3747    call wrtout(std_out,' Will try to merge this input DDB with the current one.','COLL')
3748 
3749 !  First estimate of the total number of bloks, and sto with error message if too large.
3750    write(message, '(a,i5)' ) ' Current number of bloks =',nblok
3751    call wrtout(std_out,message,'COLL')
3752    write(message, '(a,i5,a)' )' Will read ',ddb_hdr8%nblok,' blocks from the input DDB '
3753    call wrtout(std_out,message,'COLL')
3754    nblokt=nblok+ddb_hdr8%nblok
3755    if(nblokt>mblok)then
3756      write(message, '(a,i0,3a,i0,a)' )&
3757      'The expected number of blocks',nblokt,' is larger than',ch10,&
3758      'the maximum number of blocks',mblok,'.'
3759      MSG_ERROR(message)
3760    end if
3761 
3762 !  Read the bloks from the temporary database, and close it.
3763 !  Also setup the merging indicator
3764    do iblok=nblok+1,nblokt
3765      call read_blok8(ddb,iblok,ddb_hdr8%nband(1),mpert,msize,ddb_hdr8%nkpt,ddbun)
3766      mgblok(iblok)=0
3767    end do
3768    close(ddbun)
3769 
3770    nblok=nblokt
3771    write(message, '(a,i5)' ) ' Now, current number of bloks =',nblok
3772    call wrtout(std_out,message,'COLL')
3773 
3774    ! In certain cases, the different DDB will have different information
3775    ! on the pseudos (depending on fullinit)
3776    ! Here, we copy the information of the last DDB file,
3777    ! only to make the tests pass...
3778    ddb_hdr%psps%indlmn(:,:,:) = ddb_hdr8%psps%indlmn(:,:,:)
3779    ddb_hdr%psps%pspso(:) = ddb_hdr8%psps%pspso(:)
3780    ddb_hdr%psps%ekb(:,:) = ddb_hdr8%psps%ekb(:,:)
3781 
3782    call ddb_hdr_free(ddb_hdr8)
3783  end do
3784 
3785  call wrtout(std_out,' All DDBs have been read ','COLL')
3786 
3787 !*********************************************************
3788 
3789 !Check the equality of blocks, and eventually merge them
3790 
3791  if (nblok>=1) then
3792    call wrtout(std_out,' check the equality of blocks, and eventually merge ','COLL')
3793    do iblok2=2,nblok
3794      do iblok1=1,iblok2-1
3795        tmerge=0
3796 
3797 !      Check the block type identity
3798        if(ddb%typ(iblok1)==ddb%typ(iblok2))then
3799 
3800 !        Check the wavevector identities
3801          tmerge=1
3802          if(ddb%typ(iblok1)==1.or.ddb%typ(iblok1)==2)then
3803            nq=1
3804          else if(ddb%typ(iblok1)==3)then
3805 !          Note : do not merge permutation related elements ....
3806            nq=3
3807          else if(ddb%typ(iblok1)==4 .or. ddb%typ(iblok1)==0)then
3808            nq=0
3809          end if
3810          if(nq/=0)then
3811            do ii=1,nq
3812              diff= ddb%qpt(1+3*(ii-1),iblok1)/ddb%nrm(ii,iblok1) - ddb%qpt(1+3*(ii-1),iblok2)/ddb%nrm(ii,iblok2)
3813              if (abs(diff) > qtol) tmerge=0
3814              diff=ddb%qpt(2+3*(ii-1),iblok1)/ddb%nrm(ii,iblok1) - ddb%qpt(2+3*(ii-1),iblok2)/ddb%nrm(ii,iblok2)
3815              if (abs(diff) > qtol) tmerge=0
3816              diff=ddb%qpt(3+3*(ii-1),iblok1)/ddb%nrm(ii,iblok1) - ddb%qpt(3+3*(ii-1),iblok2)/ddb%nrm(ii,iblok2)
3817              if (abs(diff) > qtol) tmerge=0
3818            end do ! ii
3819          end if
3820 
3821 !        Now merges,
3822          if (tmerge == 1) then
3823            write(message, '(a,i5,a,i5)' )' merge block #',iblok2,' to block #',iblok1
3824            call wrtout(std_out,message,'COLL')
3825            mgblok(iblok2)=1
3826            do ii=1,msize
3827              if(ddb%flg(ii,iblok2)==1)then
3828                ddb%flg(ii,iblok1)=1
3829                ddb%val(1,ii,iblok1)=ddb%val(1,ii,iblok2)
3830                ddb%val(2,ii,iblok1)=ddb%val(2,ii,iblok2)
3831              end if
3832            end do
3833          end if
3834 
3835        end if
3836      end do
3837    end do
3838 
3839 !  Count the final number of bloks
3840    tmerge=0
3841    do ii=1,nblok
3842      if(mgblok(ii)==1)tmerge=tmerge+1
3843    end do
3844    nblok=nblok-tmerge
3845 
3846 !  Summarize the merging phase
3847    write(message, '(i6,a,i6,a)' )tmerge,' blocks are merged; the new DDB will have ',nblok,' blocks.'
3848    call wrtout(std_out,message,'COLL')
3849  end if !  End the condition on existence of more than one blok in current DDB
3850 
3851 !**********************************************************************
3852 
3853  write(message, '(a,a)' )' open the output database, write the',' preliminary information '
3854  call wrtout(std_out,message,'COLL')
3855 
3856  ddb_hdr%dscrpt = trim(dscrpt)
3857  ddb_hdr%nblok = nblok
3858  ddb_hdr%mblktyp = mblktyp
3859 
3860  call ddb_hdr_open_write(ddb_hdr, filnam(1), ddbun, fullinit=1)
3861 
3862  if(nddb>1)then
3863 
3864 !  Write the whole database
3865    call wrtout(std_out,' write the DDB ','COLL')
3866    choice=2
3867    do iblok=1,nblok+tmerge
3868      if(mgblok(iblok)==0)then
3869        write(std_out,'(a,i4)' ) ' Write bloc number',iblok
3870        call ddb_write_blok(ddb,iblok,choice,ddb_hdr%nband(1),mpert,msize,ddb_hdr%nkpt,ddbun)
3871      else
3872        write(message, '(a,i4,a)' )&
3873 &       ' Bloc number',iblok,' was merged, so do not write it'
3874        call wrtout(std_out,message,'COLL')
3875      end if
3876    end do
3877 
3878 !  Also write summary of bloks at the end
3879    write(ddbun, '(/,a)' )' List of bloks and their characteristics '
3880    choice=3
3881    do iblok=1,nblok+tmerge
3882      if(mgblok(iblok)==0)then
3883        call ddb_write_blok(ddb,iblok,choice,ddb_hdr%nband(1),mpert,msize,ddb_hdr%nkpt,ddbun)
3884      end if
3885    end do
3886 
3887  end if
3888 
3889  close (ddbun)
3890 
3891 !*********************************************************************
3892 
3893 !Deallocate arrays
3894 
3895  ABI_DEALLOCATE(mgblok)
3896 
3897  call ddb_hdr_free(ddb_hdr)
3898  call ddb_free(ddb)
3899 
3900 end subroutine mblktyp1

m_ddb/mblktyp5 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 mblktyp5

FUNCTION

 This routine merges the derivative databases of type 5:
 second-order eigenvalue derivatives
   why is this separate from mblktyp1? Should be merged at some point for consistency

NOTES

 The heading of the database is read, then the heading
 of the temporary database to be added is read,
 the code check their compatibility, and create a new
 database that mixes the old and the temporary ones.

TODO

  This routine is deprecated and will be removed

INPUTS

     chkopt=option for consistency checks between DDB files
     codename=MRGDDB
     ddbun=define input and output unit numbers
     dscrpt=description of the output file
     filnam=name of input or output file
     mddb=maximum number of databases (cannot be made dynamic)
     nddb=number of input DDBs

OUTPUT

     msym=maximum number of symmetry elements in space group
     Merge the file

PARENTS

      mrgddb

CHILDREN

      ddb_free,ddb_hdr_compare,ddb_hdr_free,ddb_hdr_open_read
      ddb_hdr_open_write,ddb_malloc,ddb_write_blok,read_blok8,wrtout

SOURCE

3944 subroutine mblktyp5 (chkopt,ddbun,dscrpt,filnam,mddb,msym,nddb,vrsddb)
3945 
3946 
3947 !This section has been created automatically by the script Abilint (TD).
3948 !Do not modify the following lines by hand.
3949 #undef ABI_FUNC
3950 #define ABI_FUNC 'mblktyp5'
3951 !End of the abilint section
3952 
3953  implicit none
3954 
3955 !Arguments -------------------------------
3956 !scalars
3957  integer,intent(in) :: ddbun,mddb,nddb,vrsddb
3958  integer,intent(out) :: msym
3959  character(len=fnlen),intent(in) :: dscrpt,filnam(mddb+1)
3960 
3961 !Local variables -------------------------
3962 !scalars
3963 !Define input and output unit numbers:
3964  integer,parameter :: ddbuntmp=3
3965  integer :: chkopt,choice,dimekb,iblok,iblok1,iblok2
3966  integer :: iddb,ii,lmnmax,matom
3967  integer :: mband,mblktyp,mblok,mkpt,mpert,msize,mtypat
3968  integer :: nblok,nblokt
3969  integer :: temp,tmerge,usepaw
3970  integer,allocatable :: mgblok(:)
3971  real(dp),parameter :: qtol=2.0d-8
3972  real(dp) :: diff
3973  type(ddb_type) :: ddb
3974  type(ddb_hdr_type) :: ddb_hdr, ddb_hdr8
3975 !arrays
3976  real(dp),allocatable :: blkval2(:,:,:,:),kpnt(:,:,:)
3977  character(len=500) :: message
3978 
3979 ! *********************************************************************
3980 
3981  ! Make sure there is more than one ddb to be read
3982  if(nddb==1)then
3983    write(message, '(a,a,a,a,a)' )&
3984 &   'The initialisation mode of MRGDDB, that uses nddb=1,',&
3985 &   'has been disabled in version 2.1 of ABINIT.',&
3986 &   'Action: you should use DDBs that include the symmetry',&
3987 &   'information (and that can be used and merged without',&
3988 &   'initialisation), or you should use ABINITv2.0.'
3989    MSG_ERROR(message)
3990  end if
3991 
3992 !Evaluate the maximal dimensions of arrays
3993  dimekb=0 ; matom=0 ; mband=0  ; mblok=0 ; mkpt=0
3994  msize=0  ; mtypat=0 ; lmnmax=0 ; usepaw=0 ; mblktyp = 1
3995  msym=192
3996 
3997  do iddb=1,nddb
3998    call ddb_hdr_open_read(ddb_hdr, filnam(iddb+1), ddbun, vrsddb, dimonly=1)
3999 
4000    mblok=mblok+ddb_hdr%nblok
4001    mblktyp=max(mblktyp,ddb_hdr%mblktyp)
4002    matom=max(matom,ddb_hdr%matom)
4003    mkpt=max(mkpt,ddb_hdr%mkpt)
4004    mtypat=max(mtypat,ddb_hdr%mtypat)
4005    msym=max(msym,ddb_hdr%msym)
4006    mband=max(mband,ddb_hdr%mband)
4007    dimekb=max(dimekb,ddb_hdr%psps%dimekb)
4008    lmnmax=max(lmnmax,ddb_hdr%psps%lmnmax)
4009    usepaw=max(usepaw,ddb_hdr%usepaw)
4010 
4011 
4012    call ddb_hdr_free(ddb_hdr)
4013 
4014  end do
4015 
4016  mpert=matom+6
4017  msize=3*mpert*3*mpert
4018  if(mblktyp==3)msize=msize*3*mpert
4019 
4020 !write(std_out,*),'msize',msize,'mpert',mpert,'mblktyp',mblktyp
4021  call ddb_malloc(ddb,msize,mblok,matom,mtypat)
4022 
4023 !Allocate arrays
4024  ABI_ALLOCATE(mgblok,(mblok))
4025 
4026 
4027 !**********************************************************************
4028 
4029 !Read the first database
4030 
4031  write(std_out,*)' read the input derivative database information'
4032  call ddb_hdr_open_read(ddb_hdr, filnam(2), ddbun, vrsddb, &
4033 & matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,&
4034 & msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw)
4035 
4036  ABI_ALLOCATE(blkval2,(2,msize,ddb_hdr%nband(1),mkpt))
4037  ABI_ALLOCATE(kpnt,(3,mkpt,mblok))
4038 
4039  nblok = ddb_hdr%nblok
4040 
4041  if(nblok>=1)then
4042 !  Read the blocks from the input database.
4043    write(message, '(a,i5,a)' ) ' read ',nblok,' blocks from the input DDB '
4044    call wrtout(std_out,message,'COLL')
4045    choice=1
4046    do iblok=1,nblok
4047      call read_blok8(ddb,iblok,ddb_hdr%nband(1),mpert,&
4048 &     msize,ddb_hdr%nkpt,ddbun,blkval2(1,1,1,1),kpnt(1,1,iblok))
4049 !    Setup merged indicator
4050      mgblok(iblok)=1
4051    end do
4052  else
4053    call wrtout(std_out,' No bloks in the first ddb ','COLL')
4054  end if
4055 !Close the first ddb
4056  close(ddbun)
4057 
4058 !*********************************************
4059 
4060 !In case of merging of DDBs, iterate the reading
4061  do iddb=2,nddb
4062 
4063 !  Open the corresponding input DDB,
4064 !  and read the database file information
4065    write(message, '(a,a,i6)' )ch10,&
4066 &   ' read the input derivative database number',iddb
4067    call wrtout(std_out,message,'COLL')
4068 
4069    call ddb_hdr_open_read(ddb_hdr8, filnam(iddb+1), ddbun, vrsddb, &
4070 &   matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,&
4071 &   msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw)
4072 
4073    if (chkopt==1)then
4074 !    Compare the current DDB and input DDB information.
4075 !    In case of an inconsistency, halt the execution.
4076      write(message, '(a)' )' compare the current and input DDB information'
4077      call wrtout(std_out,message,'COLL')
4078 
4079      call ddb_hdr_compare(ddb_hdr, ddb_hdr8)
4080 
4081    else if(chkopt==0)then
4082 !    No comparison between the current DDB and input DDB information.
4083      write(message, '(a)' )' no comparison between the current and input DDB information'
4084      call wrtout(std_out,message,'COLL')
4085      write(message, '(a,a,a)' )&
4086 &     'No comparison/check is performed for the current and input DDB information ',&
4087 &     'because argument --nostrict was passed to the command line. ',&
4088 &     'Use at your own risk !'
4089      MSG_COMMENT(message)
4090    end if
4091 
4092    call wrtout(std_out,' Will try to merge this input DDB with the current one.','COLL')
4093 
4094 !  First estimate of the total number of bloks, and error
4095 !  message if too large
4096    write(message, '(a,i5)' ) ' Current number of bloks =',nblok
4097    call wrtout(std_out,message,'COLL')
4098    write(message, '(a,i5,a)' )' Will read ',ddb_hdr8%nblok,' blocks from the input DDB '
4099    call wrtout(std_out,message,'COLL')
4100    nblokt=nblok+ddb_hdr8%nblok
4101    if(nblokt>mblok)then
4102      write(message, '(a,i0,a,a,a,i0,a)' )&
4103 &     'The expected number of blocks',nblokt,' is larger than',ch10,&
4104 &     'the maximum number of blocks',mblok,'.'
4105      MSG_ERROR(message)
4106    end if
4107 
4108 !  Read the bloks from the temporary database, and close it.
4109 !  Also setup the merging indicator
4110    choice=1
4111    do iblok=nblok+1,nblokt
4112      call read_blok8(ddb,iblok,ddb_hdr8%nband(1),mpert,&
4113 &     msize,ddb_hdr8%nkpt,ddbun,blkval2(1,1,1,1),kpnt(1,1,iblok))
4114      mgblok(iblok)=1
4115    end do
4116    close(ddbun)
4117 
4118    nblok=nblokt
4119    write(message, '(a,i5)' ) ' Now, current number of bloks =',nblok
4120    call wrtout(std_out,message,'COLL')
4121 
4122    ! In certain cases, the different DDB will have different information
4123    ! on the pseudos (depending on fullinit)
4124    ! Here, we copy the information of the last DDB file,
4125    ! only to make the tests pass...
4126    ddb_hdr%psps%indlmn(:,:,:) = ddb_hdr8%psps%indlmn(:,:,:)
4127    ddb_hdr%psps%pspso(:) = ddb_hdr8%psps%pspso(:)
4128    ddb_hdr%psps%ekb(:,:) = ddb_hdr8%psps%ekb(:,:)
4129 
4130    call ddb_hdr_free(ddb_hdr8)
4131 
4132  end do
4133 
4134  call wrtout(std_out,' All DDBs have been read ','COLL')
4135 
4136 !*********************************************************
4137 
4138 !Check the equality of blocks, and eventually merge them
4139 
4140  if(nblok>=1)then
4141    call wrtout(std_out,' check the equality of blocks, and eventually merge ','COLL')
4142    do iblok2=2,nblok
4143      do iblok1=1,iblok2-1
4144 !      Check the block type identity
4145        if(ddb%typ(iblok1)==ddb%typ(iblok2))then
4146 !        Check the wavevector identities
4147          diff=abs(ddb%qpt(1,iblok1)-ddb%qpt(1,iblok2))
4148          diff=diff+abs(ddb%qpt(2,iblok1)-ddb%qpt(2,iblok2))
4149          diff=diff+abs(ddb%qpt(3,iblok1)-ddb%qpt(3,iblok2))
4150          if(abs(diff)<qtol)mgblok(iblok2)=0
4151        end if
4152      end do
4153    end do
4154 
4155 !  Count the final number of bloks
4156    tmerge=0
4157    do ii=1,nblok
4158      if(mgblok(ii)==1)tmerge=tmerge+1
4159    end do
4160    temp = nblok-tmerge
4161    nblok=tmerge
4162 
4163 !  Summarize the merging phase
4164    write(message, '(i6,a,i6,a)' )&
4165 &   temp,' blocks are merged; the new DDB will have ',nblok,' blocks.'
4166    call wrtout(std_out,message,'COLL')
4167  end if
4168 
4169 !**********************************************************************
4170 
4171  write(message, '(a,a)' )' open the output database, write the',' preliminary information '
4172  call wrtout(std_out,message,'COLL')
4173 
4174  ddb_hdr%dscrpt = trim(dscrpt)
4175  ddb_hdr%nblok = nblok !nblokt
4176  ddb_hdr%mblktyp = mblktyp
4177 
4178  call ddb_hdr_open_write(ddb_hdr, filnam(1), ddbun, fullinit=1)
4179 
4180  if(nddb>1)then
4181 
4182 !  Write the whole database
4183    call wrtout(std_out,' write the DDB ','COLL')
4184    ii = 1 !unit indicator of what will be merged
4185 !  Create a temporary file to decrease memory need.
4186    do iddb=1,nddb
4187      call ddb_hdr_open_read(ddb_hdr8, filnam(iddb+1), ddbuntmp, vrsddb)
4188 
4189      do iblok=1,ddb_hdr8%nblok
4190        if(mgblok(ii)==1) then
4191          call read_blok8(ddb,ii,ddb_hdr8%nband(1),mpert,&
4192 &         msize,ddb_hdr8%nkpt,ddbuntmp,blkval2(:,:,:,:),kpnt(:,:,ii))
4193          choice=2
4194          call ddb_write_blok(ddb,ii,choice,ddb_hdr%nband(1),mpert,&
4195 &         msize,ddb_hdr8%nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,ii))
4196        else
4197          write(message, '(a,i4,a,i4,a)' )&
4198 &         ' Bloc number',iblok,' of DDB ',iddb,&
4199 &         ' was merged, so do not write it'
4200          call wrtout(std_out,message,'COLL')
4201        end if
4202        ii = ii+1
4203      end do
4204      close(ddbuntmp)
4205      call ddb_hdr_free(ddb_hdr8)
4206    end do !iddb=1,nddb
4207 
4208 !  Also write summary of bloks at the end
4209    write(ddbun, '(/,a)' )' List of bloks and their characteristics '
4210    choice=3
4211    do iblok=1,nblokt
4212      if(mgblok(iblok)==1)then
4213        call ddb_write_blok(ddb,iblok,choice,ddb_hdr%nband(1),mpert,&
4214 &       msize,ddb_hdr%nkpt,ddbun)
4215      end if
4216    end do
4217 
4218  end if
4219 
4220  close (ddbun)
4221 
4222 !*********************************************************************
4223 
4224  call ddb_hdr_free(ddb_hdr)
4225 
4226 !Deallocate arrays
4227  ABI_DEALLOCATE(mgblok)
4228  ABI_DEALLOCATE(blkval2)
4229  ABI_DEALLOCATE(kpnt)
4230 
4231  call ddb_free(ddb)
4232 
4233 end subroutine mblktyp5

m_ddb/nlopt [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 nlopt

FUNCTION

 Output of all quantities related to third-order derivatives of the energy.
 Compute the permutations of the three perturbations, then
 write out the whole matrix of third order derivatives
 in reduced coordinates. Finally, compute the non-linear optical
 susceptibility d and the first-order change in the dielectric
 susceptibility tensor induced by an atomic displacement.

INPUTS

  blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte
   has been calculated ; 0 otherwise )
  d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
  gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
  mpert =maximum number of ipert
  natom= number of atoms
  rprimd(3,3)=dimensional primitive translations (bohr)
  ucvol=unit cell volume (bohr^3)

OUTPUT

 carflg(3,mpert,3,mpert,3,mpert)=1 if the element of d3cart has been calculated, 0 otherwise
 d3cart(2,3,mpert,3,mpert,3,mpert)=matrix of third-order energy derivatives in cartesian coordinates

PARENTS

      m_ddb,nonlinear

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1473 subroutine nlopt(blkflg,carflg,d3,d3cart,gprimd,mpert,natom,rprimd,ucvol)
1474 
1475 
1476 !This section has been created automatically by the script Abilint (TD).
1477 !Do not modify the following lines by hand.
1478 #undef ABI_FUNC
1479 #define ABI_FUNC 'nlopt'
1480 !End of the abilint section
1481 
1482  implicit none
1483 
1484 !Arguments -------------------------------
1485 !scalars
1486  integer,intent(in) :: mpert,natom
1487  real(dp),intent(in) :: ucvol
1488 !arrays
1489  integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert)
1490  integer,intent(out) :: carflg(3,mpert,3,mpert,3,mpert)
1491  real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert),gprimd(3,3),rprimd(3,3)
1492  real(dp),intent(out) :: d3cart(2,3,mpert,3,mpert,3,mpert)
1493 
1494 !Local variables -------------------------
1495 !scalars
1496  integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
1497 !arrays
1498  integer :: flg1(3),flg2(3)
1499  real(dp) :: vec1(3),vec2(3)
1500 
1501 ! *******************************************************************
1502 
1503 !Compute the permutations of the perturbations
1504 
1505  d3cart(:,:,:,:,:,:,:) = 0._dp
1506 
1507  do i1pert = 1,mpert
1508    do i2pert = 1,mpert
1509      do i3pert = 1,mpert
1510        do i1dir=1,3
1511          do i2dir=1,3
1512            do i3dir=1,3
1513 
1514 !            Check if all elements are available
1515 
1516              if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=0).and. &
1517 &             (blkflg(i1dir,i1pert,i3dir,i3pert,i2dir,i2pert)/=0).and. &
1518 &             (blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert)/=0).and. &
1519 &             (blkflg(i2dir,i2pert,i3dir,i3pert,i1dir,i1pert)/=0).and. &
1520 &             (blkflg(i3dir,i3pert,i1dir,i1pert,i2dir,i2pert)/=0).and. &
1521 &             (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=0)) then
1522 
1523                d3cart(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
1524 &               (  d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + &
1525 &               d3(:,i1dir,i1pert,i3dir,i3pert,i2dir,i2pert) + &
1526 &               d3(:,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) + &
1527 &               d3(:,i2dir,i2pert,i3dir,i3pert,i1dir,i1pert) + &
1528 &               d3(:,i3dir,i3pert,i1dir,i1pert,i2dir,i2pert) + &
1529 &               d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert))*sixth
1530 
1531              end if
1532            end do
1533          end do
1534        end do
1535      end do
1536    end do
1537  end do
1538 
1539 !Transform to cartesian coordinates
1540  carflg(:,:,:,:,:,:) = 0
1541 
1542  do i1pert = 1, mpert
1543    do i2pert = 1, mpert
1544      do i3pert = 1, mpert
1545 
1546        do i2dir = 1, 3
1547          do i3dir = 1, 3
1548 
1549            vec1(:) = d3cart(1,:,i1pert,i2dir,i2pert,i3dir,i3pert)
1550            flg1(:) = blkflg(:,i1pert,i2dir,i2pert,i3dir,i3pert)
1551            call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2)
1552            d3cart(1,:,i1pert,i2dir,i2pert,i3dir,i3pert) = vec2(:)
1553            carflg(:,i1pert,i2dir,i2pert,i3dir,i3pert) = flg2(:)
1554 
1555          end do
1556        end do
1557 
1558        do i1dir = 1, 3
1559          do i3dir = 1, 3
1560            vec1(:) = d3cart(1,i1dir,i1pert,:,i2pert,i3dir,i3pert)
1561            flg1(:) = blkflg(i1dir,i1pert,:,i2pert,i3dir,i3pert)
1562            call cart39(flg1,flg2,gprimd,i2pert,natom,rprimd,vec1,vec2)
1563            d3cart(1,i1dir,i1pert,:,i2pert,i3dir,i3pert) = vec2(:)
1564            carflg(i1dir,i1pert,:,i2pert,i3dir,i3pert) = flg2(:)
1565          end do
1566        end do
1567 
1568        do i1dir = 1, 3
1569          do i2dir = 1, 3
1570            vec1(:) = d3cart(1,i1dir,i1pert,i2dir,i2pert,:,i3pert)
1571            flg1(:) = blkflg(i1dir,i1pert,i2dir,i2pert,:,i3pert)
1572            call cart39(flg1,flg2,gprimd,i3pert,natom,rprimd,vec1,vec2)
1573            d3cart(1,i1dir,i1pert,i2dir,i2pert,:,i3pert) = vec2(:)
1574            carflg(i1dir,i1pert,i2dir,i2pert,:,i3pert) = flg2(:)
1575          end do
1576        end do
1577 
1578      end do
1579    end do
1580  end do
1581 
1582 !Compute non linear-optical coefficients d_ijk (atomic units)
1583 
1584  i1pert = natom+2
1585  d3cart(:,:,i1pert,:,i1pert,:,i1pert) = -3._dp*d3cart(:,:,i1pert,:,i1pert,:,i1pert)/(ucvol*2._dp)
1586 
1587 !Compute first-order change in the electronic dielectric
1588 !susceptibility (Bohr^-1) induced by an atomic displacement
1589  d3cart(1:2,1:3,1:natom,1:3,natom + 2,1:3,natom + 2) = &
1590 & -6._dp*d3cart(1:2,1:3,1:natom,1:3,natom + 2,1:3,natom + 2)/ucvol
1591 
1592 end subroutine nlopt

m_ddb/rdddb9 [ Functions ]

[ Top ] [ m_ddb ] [ Functions ]

NAME

 rdddb9

FUNCTION

 This routine reads the derivative database entirely,
 for use in ppddb9, and performs some checks and symmetrisation
 At the end, the whole DDB is in central memory, contained in the
 array ddb%val(2,msize,ddb%nblok).
 The information on it is contained in the four arrays
   ddb%flg(msize,ddb%nblok) : blok flag for each element
   ddb%qpt(9,ddb%nblok)  : blok wavevector (unnormalized)
   ddb%nrm(3,ddb%nblok)  : blok wavevector normalization
   ddb%typ(ddb%nblok)    : blok type

INPUTS

 atifc(natom) = atifc(ia) equals 1 if the analysis of ifc has to be done for atom ia; otherwise 0
 ddbun = unit number for DDB io
 dimekb=dimension of ekb (for the time being, only for norm- conserving psps)
 iout=unit number for output of formatted data
 filnam=name of input file
 lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
       =if useylm=0, max number of (l,n)   comp. over all type of psps
 mband=maximum number of bands
 mpert =maximum number of ipert
 msize=maximum size of data blocks
 msym =maximum number of symmetry elements in space group
 natifc = number of atoms for which the analysis of ifc is done
 natom = number of atoms
 ntypat=number of atom types
 usepaw= 0 for non paw calculation; =1 for paw calculation

OUTPUT

 acell(3)=length scales of cell (bohr)
 amu(ntypat)=mass of the atoms (atomic mass unit)
 ddb : ddb blok datatype
   contents: ddb%flg(msize,nblok)= flag of existence for each element of the DDB
             ddb%nrm(3,nblok)  : blok wavevector normalization
             ddb%qpt(9,nblok)  : blok wavevector (unnormalized)
             ddb%typ(nblok)    : blok type
             ddb%val(2,msize,nblok)= value of each complex element of the DDB
             ddb%nblok= number of bloks in the DDB
 gmet(3,3)=reciprocal space metric tensor in bohr**-2
 gprim(3,3)=dimensionless reciprocal space primitive translations
 indsym(4,msym,natom)=indirect indexing array for symmetries
 natom=number of atoms in cell
 nsym=number of space group symmetries
 rmet(3,3)=metric tensor in real space (bohr^2)
 rprim(3,3)= primitive translation vectors
 symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
 symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
 symafm(nsym)=Anti-ferromagnetic symmetries.
 tnons(3,nsym)=fractional nonsymmorphic translations
 typat(natom)=type integer for each atom in cell
 ucvol=unit cell volume in bohr**3
 xcart(3,natom)=atomic cartesian coordinates
 xred(3,natom)=fractional dimensionless atomic coordinates
 zion(ntypat)=charge on each type of atom (real number)
 znucl(ntypat)=Nuclear charge for each type of pseudopotential

PARENTS

      m_ddb

CHILDREN

      ddb_hdr_free,ddb_hdr_open_read

SOURCE

1134 subroutine rdddb9(acell,atifc,amu,ddb,ddbun,filnam,gmet,gprim,indsym,iout,&
1135 & mband,mpert,msize,msym,natifc,natom,nkpt,nsym,ntypat,&
1136 & rmet,rprim,symrec,symrel,symafm,tnons,typat,ucvol,xcart,xred,zion,znucl)
1137 
1138 
1139 !This section has been created automatically by the script Abilint (TD).
1140 !Do not modify the following lines by hand.
1141 #undef ABI_FUNC
1142 #define ABI_FUNC 'rdddb9'
1143 !End of the abilint section
1144 
1145  implicit none
1146 
1147 !Arguments ------------------------------------
1148 ! NOTE: these are used for dimensioning and then re-assigned in ioddb8.
1149 !   This is almost definitely bad practice. In particular
1150 !    it should be indsym(4,msym,natom),
1151 !   and
1152 !    the allocation allocate(kpt(3,nkpt)) is strange
1153 !scalars
1154  integer,intent(in) :: ddbun,iout,mband,mpert,msize,msym,natifc
1155  integer,intent(inout) :: natom,nkpt,nsym,ntypat
1156  real(dp),intent(out) :: ucvol
1157  character(len=*),intent(in) :: filnam
1158  type(ddb_type),intent(inout) :: ddb
1159 !arrays
1160  integer,intent(inout) :: atifc(natom)
1161  integer,intent(inout) :: indsym(4,msym,natom)
1162  integer,intent(out) :: symrec(3,3,msym),symrel(3,3,msym),symafm(msym)
1163  integer,intent(out) :: typat(natom)
1164  real(dp),intent(out) :: acell(3),amu(ntypat)
1165  real(dp),intent(out) :: gmet(3,3),gprim(3,3),rmet(3,3)
1166  real(dp),intent(out) :: rprim(3,3),tnons(3,msym),xcart(3,natom),xred(3,natom)
1167  real(dp),intent(out) :: zion(ntypat),znucl(ntypat)
1168 
1169 !Local variables -------------------------
1170 !mtyplo=maximum number of type, locally
1171 !scalars
1172  integer,parameter :: msppol=2,mtyplo=6
1173  integer :: iblok,isym
1174  integer :: nsize,timrev
1175  integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
1176  real(dp),parameter :: tolsym8=tol8
1177  character(len=500) :: message
1178  type(ddb_hdr_type) :: ddb_hdr
1179 !arrays
1180  integer :: symq(4,2,msym)
1181  integer,allocatable :: car3flg(:,:,:,:,:,:),carflg(:,:,:,:)
1182  integer,allocatable :: tmpflg(:,:,:,:,:,:),rfpert(:,:,:,:,:,:)
1183  real(dp) :: gprimd(3,3),qpt(3),rprimd(3,3)
1184  real(dp),allocatable :: d2cart(:,:,:,:,:),d3cart(:,:,:,:,:,:,:)
1185  real(dp),allocatable :: tmpval(:,:,:,:,:,:,:)
1186 
1187 ! *********************************************************************
1188 
1189  DBG_ENTER("COLL")
1190 
1191 !Open the input derivative database file and read the header
1192  call ddb_hdr_open_read(ddb_hdr, filnam, ddbun, DDB_VERSION, msym=msym, mband=mband)
1193 
1194  !nkpt = ddb_hdr%nkpt
1195  !ntypat = ddb_hdr%ntypat
1196  nsym = ddb_hdr%nsym
1197  acell = ddb_hdr%acell
1198  rprim = ddb_hdr%rprim
1199 
1200  amu(:) = ddb_hdr%amu(1:ntypat)
1201  typat(:) = ddb_hdr%typat(1:natom)
1202  zion(:) = ddb_hdr%zion(1:ntypat)
1203  znucl(:) = ddb_hdr%znucl(1:ntypat)
1204 
1205  symafm(:) = ddb_hdr%symafm(:)
1206  symrel(:,:,:) = ddb_hdr%symrel(:,:,:)
1207  tnons(:,:) = ddb_hdr%tnons(:,:)
1208 
1209  xred(:,:) = ddb_hdr%xred(:,:)
1210 
1211  call ddb_hdr_free(ddb_hdr)
1212 
1213 !Compute different matrices in real and reciprocal space, also
1214 !checks whether ucvol is positive.
1215  call mkrdim(acell,rprim,rprimd)
1216  call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
1217 
1218 !Obtain reciprocal space primitive transl g from inverse trans of r
1219 !(Unlike in abinit, gprim is used throughout ifc; should be changed, later)
1220  call matr3inv(rprim,gprim)
1221 
1222 !Generate atom positions in cartesian coordinates
1223  call xred2xcart(natom,rprimd,xcart,xred)
1224 
1225 !Transposed inversion of the symmetry matrices, for use in the reciprocal space
1226  do isym=1,nsym
1227    call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
1228  end do
1229 
1230 !SYMATM generates for all the atoms and all the symmetries, the atom
1231 !on which the referenced one is sent and also the translation bringing
1232 !back this atom to the referenced unit cell
1233  call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred)
1234 
1235 !Check the correctness of some input parameters, and perform small treatment if needed.
1236  call chkin9(atifc,natifc,natom)
1237 
1238 !Read the blocks from the input database, and close it.
1239  write(message, '(3a,i0,a)' )ch10,ch10,' rdddb9: read ',ddb%nblok,' blocks from the input DDB '
1240  call wrtout(std_out,message,'COLL')
1241 
1242  do iblok=1,ddb%nblok
1243    call read_blok8(ddb,iblok,mband,mpert,msize,nkpt,ddbun)
1244 
1245    !  Here complete the matrix by symmetrisation of the existing elements
1246    if(ddb%typ(iblok)==1 .or. ddb%typ(iblok)==2) then
1247 
1248      qpt(1)=ddb%qpt(1,iblok)/ddb%nrm(1,iblok)
1249      qpt(2)=ddb%qpt(2,iblok)/ddb%nrm(1,iblok)
1250      qpt(3)=ddb%qpt(3,iblok)/ddb%nrm(1,iblok)
1251 
1252      ! Examine the symmetries of the q wavevector
1253      call littlegroup_q(nsym,qpt,symq,symrec,symafm,timrev,prtvol=0)
1254 
1255      nsize=3*mpert*3*mpert
1256      ABI_MALLOC(tmpflg,(3,mpert,3,mpert,1,1))
1257      ABI_MALLOC(tmpval,(2,3,mpert,3,mpert,1,1))
1258 
1259      tmpflg(:,:,:,:,1,1) = reshape(ddb%flg(1:nsize,iblok), shape = (/3,mpert,3,mpert/))
1260      tmpval(1,:,:,:,:,1,1) = reshape(ddb%val(1,1:nsize,iblok), shape = (/3,mpert,3,mpert/))
1261      tmpval(2,:,:,:,:,1,1) = reshape(ddb%val(2,1:nsize,iblok), shape = (/3,mpert,3,mpert/))
1262 
1263      ! Then apply symmetry operations
1264      call d2sym3(tmpflg,tmpval,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev,1)
1265 
1266      ! Transform the dynamical matrix in cartesian coordinates
1267      ABI_MALLOC(carflg,(3,mpert,3,mpert))
1268      ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
1269 
1270      call cart29(tmpflg,tmpval,carflg,d2cart,gprimd,1,mpert,natom,1,ntypat,rprimd,typat,ucvol,zion)
1271 
1272      ddb%flg(1:nsize,iblok) = reshape(carflg,shape = (/3*mpert*3*mpert/))
1273      ddb%val(1,1:nsize,iblok) = reshape(d2cart(1,:,:,:,:), shape = (/3*mpert*3*mpert/))
1274      ddb%val(2,1:nsize,iblok) = reshape(d2cart(2,:,:,:,:), shape = (/3*mpert*3*mpert/))
1275 
1276      ABI_FREE(carflg)
1277      ABI_FREE(d2cart)
1278      ABI_FREE(tmpflg)
1279      ABI_FREE(tmpval)
1280 
1281    else if (ddb%typ(iblok) == 3) then
1282 
1283      nsize=3*mpert*3*mpert*3*mpert
1284      ABI_MALLOC(tmpflg,(3,mpert,3,mpert,3,mpert))
1285      ABI_MALLOC(tmpval,(2,3,mpert,3,mpert,3,mpert))
1286      ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert))
1287 
1288      tmpflg(:,:,:,:,:,:) = reshape(ddb%flg(1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/))
1289      tmpval(1,:,:,:,:,:,:) = reshape(ddb%val(1,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/))
1290      tmpval(2,:,:,:,:,:,:) = reshape(ddb%val(2,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/))
1291 
1292 !    Set the elements that are zero by symmetry for raman and
1293 !    non-linear optical susceptibility tensors
1294      rfpert = 0
1295      rfpert(:,natom+2,:,natom+2,:,natom+2) = 1
1296      rfpert(:,1:natom,:,natom+2,:,natom+2) = 1
1297      rfpert(:,natom+2,:,1:natom,:,natom+2) = 1
1298      rfpert(:,natom+2,:,natom+2,:,1:natom) = 1
1299      call sytens(indsym,mpert,natom,nsym,rfpert,symrec,symrel)
1300      do i1pert = 1,mpert
1301        do i2pert = 1,mpert
1302          do i3pert = 1,mpert
1303            do i1dir=1,3
1304              do i2dir=1,3
1305                do i3dir=1,3
1306                  if ((rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) .and. &
1307 &                  (tmpflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)) then
1308                    tmpval(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
1309                    tmpflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=1
1310                  end if
1311                end do
1312              end do
1313            end do
1314          end do
1315        end do
1316      end do
1317 
1318      call d3sym(tmpflg,tmpval,indsym,mpert,natom,nsym,symrec,symrel)
1319 
1320      ABI_MALLOC(d3cart,(2,3,mpert,3,mpert,3,mpert))
1321      ABI_MALLOC(car3flg,(3,mpert,3,mpert,3,mpert))
1322 
1323      call nlopt(tmpflg,car3flg,tmpval,d3cart,gprimd,mpert,natom,rprimd,ucvol)
1324 
1325      ddb%flg(1:nsize,iblok) = reshape(car3flg, shape = (/3*mpert*3*mpert*3*mpert/))
1326      ddb%val(1,1:nsize,iblok) = reshape(d3cart(1,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/))
1327      ddb%val(2,1:nsize,iblok) = reshape(d3cart(2,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/))
1328 
1329      ABI_FREE(d3cart)
1330      ABI_FREE(car3flg)
1331      ABI_FREE(tmpflg)
1332      ABI_FREE(tmpval)
1333      ABI_FREE(rfpert)
1334    end if
1335  end do ! iblok
1336 
1337  close(ddbun)
1338 
1339  write(message,'(a)' )' Now the whole DDB is in central memory '
1340  call wrtout(std_out,message,'COLL')
1341  call wrtout(iout,message,'COLL')
1342 
1343  DBG_EXIT("COLL")
1344 
1345 end subroutine rdddb9