TABLE OF CONTENTS


ABINIT/m_anharmonics_terms [ Functions ]

[ Top ] [ Functions ]

NAME

 m_anharmonics_term

FUNCTION

 Module with datatype and tools for the anharmonics terms

COPYRIGHT

 Copyright (C) 2010-2018 ABINIT group (AM)
 This file is distributed under the terms of the
 GNU General Public Licence, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_anharmonics_terms
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_polynomial_coeff
32  use m_ifc, only : ifc_type,ifc_free
33 
34  implicit none
35 
36  public :: anharmonics_terms_init
37  public :: anharmonics_terms_free
38  public :: anharmonics_terms_freeCoeffs
39  public :: anharmonics_terms_evaluateElastic
40  public :: anharmonics_terms_evaluateIFCStrainCoupling
41  public :: anharmonics_terms_setCoeffs
42  public :: anharmonics_terms_setElastic3rd
43  public :: anharmonics_terms_setElastic4th
44  public :: anharmonics_terms_setElasticDispCoupling
45  public :: anharmonics_terms_setStrainPhononCoupling

defs_abitypes/anharmonics_terms_type [ Types ]

[ Top ] [ defs_abitypes ] [ Types ]

NAME

 anharmonics_terms_type

FUNCTION

 datatype for a effective potential constructed.

SOURCE

58  type, public :: anharmonics_terms_type
59 
60    integer :: ncoeff = 0
61 !    nterm store the number of coefficients
62 
63    logical ::  has_elastic3rd
64 !   Flag to know if the 3rd derivatives with respect to strain is present
65 
66    logical ::  has_elastic4th
67 !   Flag to know if the 3rd derivatives with respect to strain is present
68 
69    logical ::  has_strain_coupling
70 !   Flag to know if the 3rd derivatives with respect to strain and 2 atom disp is present
71 
72    logical ::  has_elastic_displ
73 !   Flag to know if the 3rd derivatives with respect to 2 strain and 3 atom disp is present
74 
75    logical :: bounded
76 !   True : the model is bounded
77 
78    type(polynomial_coeff_type),dimension(:),allocatable :: coefficients
79 !    array with all the coefficients from  polynomial coefficients
80 
81    real(dp) :: elastic3rd(6,6,6)
82 !    elastic_constant(6,6,6)
83 !    Elastic tensor Hartree
84 
85    real(dp) :: elastic4th(6,6,6,6)
86 !    elastic_constant(6,6,6)
87 !    Elastic tensor Hartree
88 
89    real(dp), allocatable :: elastic_displacement(:,:,:,:)
90 !    elastic_displacement(6,6,3,natom)
91 !    internal strain tensor
92 
93    type(ifc_type),dimension(:),allocatable :: phonon_strain
94 !   Array of ifc with phonon_strain coupling for each strain
95 
96  end type anharmonics_terms_type

m_anharmonics_terms/anharmonics_terms_evaluateIFCStrainCoupling [ Functions ]

[ Top ] [ m_anharmonics_terms ] [ Functions ]

NAME

  anharmonics_terms_evaluateIFCStrainCoupling

FUNCTION

  This fonction compute the harmonic part of the energy
  of the supercell in the eff_pot

INPUTS

  strain_phonon(6)<type(ifc_type) = strain-phonon coupling
  disp(3,natom_sc) = atomics displacement between configuration and the reference
  natom = number of atoms in the supercell
  natom_uc = number of atoms in the unit cell
  sc_size(3) = size of the supercell
  cells(ncell) = number of the cells into the supercell (1,2,3,4,5)
  ncell  = total number of cell to treat
  index_cells(3,ncell) = indexes of the cells into  supercell (-1 -1 -1 ,...,1 1 1)
  comm=MPI communicator

OUTPUT

   energy = contribution of the ifc to the energy
   fcart(3,natom) = contribution of the ifc to the forces
   strten(6) = contribution to the stress tensor

 PARENT
   effective_potential_evaluate

PARENTS

      m_effective_potential

CHILDREN

      getpbcindexes_supercell,xmpi_sum

SOURCE

860 subroutine anharmonics_terms_evaluateIFCStrainCoupling(phonon_strain,disp,energy,fcart,natom,natom_uc,&
861 &                                                      sc_size,strain,strten,cells,ncell,&
862 &                                                      index_cells,comm)
863 
864 
865 !This section has been created automatically by the script Abilint (TD).
866 !Do not modify the following lines by hand.
867 #undef ABI_FUNC
868 #define ABI_FUNC 'anharmonics_terms_evaluateIFCStrainCoupling'
869 !End of the abilint section
870 
871  implicit none
872 
873 !Arguments -------------------------------
874 ! scalars
875   real(dp),intent(out) :: energy
876   integer,intent(in) :: natom,natom_uc,ncell
877   integer,intent(in) :: comm
878 ! array
879   integer,intent(in) ::   cells(ncell),index_cells(ncell,3)
880   integer,intent(in) :: sc_size(3)
881   type(ifc_type),intent(in) :: phonon_strain(6)
882   real(dp),intent(in) :: disp(3,natom)
883   real(dp),intent(out) :: fcart(3,natom)
884   real(dp),intent(out) :: strten(6)
885   real(dp),intent(in) :: strain(6)
886 !Local variables-------------------------------
887 ! scalar
888   integer :: alpha
889   integer :: i1,i2,i3,ia,ib,icell,ii
890   integer :: irpt,jj,kk,ll,mu,nu
891   integer :: ierr
892   real(dp):: ifc
893 ! array
894   integer :: cell_atom2(3)
895   character(500) :: msg
896 
897 ! *************************************************************************
898 
899   if (any(sc_size <= 0)) then
900     write(msg,'(a,a)')' sc_size can not be inferior or equal to zero'
901     MSG_ERROR(msg)
902   end if
903 
904 ! Initialisation of variables
905   energy   = zero
906   fcart(:,:) = zero
907   strten(:) = zero
908 
909   do icell = 1,ncell
910     ii = (cells(icell)-1)*natom_uc
911     i1=index_cells(icell,1); i2=index_cells(icell,2); i3=index_cells(icell,3)
912     do alpha=1,6
913       do irpt = 1,phonon_strain(alpha)%nrpt
914 !       get the cell of atom2  (0 0 0, 0 0 1...)
915         cell_atom2(1) =  i1 + phonon_strain(alpha)%cell(1,irpt)
916         cell_atom2(2) =  i2 + phonon_strain(alpha)%cell(2,irpt)
917         cell_atom2(3) =  i3 + phonon_strain(alpha)%cell(3,irpt)
918         call getPBCIndexes_supercell(cell_atom2(1:3),sc_size(1:3))
919 !       index of the second atom in the displacement array
920         jj = (cell_atom2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
921 &            (cell_atom2(2)-1)*sc_size(3)*natom_uc+&
922 &            (cell_atom2(3)-1)*natom_uc
923         do ib = 1, natom_uc
924           ll = jj + ib
925           do nu=1,3
926             do ia = 1, natom_uc
927               kk = ii + ia
928               do mu=1,3
929                 ifc = phonon_strain(alpha)%atmfrc(mu,ia,nu,ib,irpt)
930 !               accumule energy
931                 energy =  energy + sixth*strain(alpha)*disp(mu,kk)*disp(nu,ll)*ifc
932 !               accumule forces
933                 fcart(mu,kk) = fcart(mu,kk) + half*strain(alpha)*disp(nu,ll)*ifc
934 !               accumule stresses
935                 strten(alpha) = strten(alpha) + half*disp(mu,kk)*disp(nu,ll)*ifc
936               end do
937             end do
938           end do
939         end do
940       end do
941     end do
942   end do
943 
944 ! MPI_SUM
945   call xmpi_sum(energy, comm, ierr)
946   call xmpi_sum(fcart , comm, ierr)
947   call xmpi_sum(strten, comm, ierr)
948 
949 end subroutine anharmonics_terms_evaluateIFCStrainCoupling

m_anharmonics_terms/anharmonics_terms_init [ Functions ]

[ Top ] [ m_anharmonics_terms ] [ Functions ]

NAME

 anharmonics_terms_init

FUNCTION

 Initialize anharmonics_terms datatype

INPUTS

 natom  = number of atoms in primitive cell
 ncoeff = number of coefficient for the fited polynome
 bounded = optional, flag to now if the model in bounded
 elastic3rd(6,6,6) = optional,3rd order of the elastic constants
 elastic4th(6,6,6,6) = optional,4st order of the elastic constants
 elastic_displacement(6,6,3,natom) = optional,elastic constant - force coupling
 phonon_strain<type(ifc_type)>(6) = optional,phonon strain couling
 coeffs<type(polynomial_coeff_type)>(ncoeff) = optional,datatype with polynomial coefficients

OUTPUT

 anharmonics_terms<type(anharmonics_terms_type)> = anharmonics_terms datatype to be initialized

PARENTS

      m_effective_potential

CHILDREN

      getpbcindexes_supercell,xmpi_sum

SOURCE

131 subroutine anharmonics_terms_init(anharmonics_terms,natom,ncoeff,&
132 &                                 bounded,elastic3rd,elastic4th,elastic_displacement,&
133 &                                 phonon_strain,coeffs)
134 
135 
136 !This section has been created automatically by the script Abilint (TD).
137 !Do not modify the following lines by hand.
138 #undef ABI_FUNC
139 #define ABI_FUNC 'anharmonics_terms_init'
140 !End of the abilint section
141 
142  implicit none
143 
144 !Arguments ------------------------------------
145 !scalars
146  integer, intent(in) :: natom,ncoeff
147  type(anharmonics_terms_type), intent(out) :: anharmonics_terms
148  real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
149  real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
150  type(polynomial_coeff_type),optional :: coeffs(ncoeff)
151  type(ifc_type),optional,intent(in) :: phonon_strain(6)
152  logical,optional,intent(in) :: bounded
153 !arrays
154 !Local variables-------------------------------
155 !scalar
156 !arrays
157  character(len=500) :: msg
158 
159 ! *************************************************************************
160 
161  call anharmonics_terms_free(anharmonics_terms)
162 
163 ! Check the number of atoms
164  if (natom < 1) then
165    write(msg, '(a,a,a,i10,a)' )&
166 &   'The cell must have at least one atom.',ch10,&
167 &   'The number of atom is  ',natom,'.'
168    MSG_BUG(msg)
169  end if
170 
171 !Allocation of phonon strain coupling array (3rd order)
172  if(present(phonon_strain)) then
173    call anharmonics_terms_setStrainPhononCoupling(anharmonics_terms,natom,phonon_strain)
174  end if
175 
176 !Set the 3rd order elastic tensor
177  anharmonics_terms%elastic3rd = zero
178  if(present(elastic3rd))then
179    call anharmonics_terms_setElastic3rd(anharmonics_terms,elastic3rd)
180  end if
181 
182 !Set the 3rd order elastic tensor
183  anharmonics_terms%elastic4th = zero
184  if(present(elastic4th))then
185    call anharmonics_terms_setElastic4th(anharmonics_terms,elastic4th)
186  end if
187 
188 !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement
189  if(present(elastic_displacement))then
190    call anharmonics_terms_setElasticDispCoupling(anharmonics_terms,natom,elastic_displacement)
191  end if
192 
193  anharmonics_terms%ncoeff = 0
194 
195 !Allocation of the coefficient
196   if(present(coeffs))then
197    if(ncoeff /= size(coeffs))then
198      write(msg, '(a)' )&
199 &        ' ncoeff has not the same size than coeffs array, '
200      MSG_BUG(msg)
201    end if
202    call anharmonics_terms_setCoeffs(coeffs,anharmonics_terms,ncoeff)
203  end if
204 
205 !Set the flag bounded
206  if(present(bounded))then
207    anharmonics_terms%bounded = bounded
208  else
209    anharmonics_terms%bounded = .FALSE.
210  end if
211 
212 end subroutine anharmonics_terms_init

m_effective_potential/anharmonics_terms_evaluateElastic [ Functions ]

[ Top ] [ m_effective_potential ] [ Functions ]

NAME

  anharmonics_terms_evaluateElastic

FUNCTION

 Compute the energy, stresses and forces related to the application of strain

INPUTS

 disp(3,natom_sc) = atomics displacement between configuration and the reference
 natom = number of atom in the supercell
 natom_uc = number of atom in the unit cell
 ncell  = number of cell
 strain(6) =  strain to apply
 elastic3rd(6,6,6) = 3 order derivatives with respect to to 3 strain
 elastic4th(6,6,66,) = 4 order derivatives with respect to to 4 strain
 elastic_displacement(6,6,3,natom) = 3 order derivatives with respect to 2 strain and 1 Atom disp

OUTPUT

   energy = contribution of the ifc to the energy
   fcart(3,natom) = contribution of the ifc to the forces
   strten(6) = contribution to the stress tensor

PARENTS

      m_effective_potential

CHILDREN

      getpbcindexes_supercell,xmpi_sum

SOURCE

733 subroutine anharmonics_terms_evaluateElastic(disp,energy,fcart,natom,natom_uc,ncell,strten,strain,&
734 &                                            elastic3rd,elastic4th,elastic_displacement)
735 
736 
737 !This section has been created automatically by the script Abilint (TD).
738 !Do not modify the following lines by hand.
739 #undef ABI_FUNC
740 #define ABI_FUNC 'anharmonics_terms_evaluateElastic'
741 !End of the abilint section
742 
743  real(dp),intent(out):: energy
744  integer, intent(in) :: natom,natom_uc,ncell
745 ! array
746  real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
747  real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
748  real(dp),intent(out):: strten(6)
749  real(dp),intent(out):: fcart(3,natom)
750  real(dp),intent(in) :: disp(3,natom)
751  real(dp),intent(in) :: strain(6)
752 
753  !Local variables-------------------------------
754 ! scalar
755  integer :: ia,ii,mu,alpha,beta,gamma,delta,d1,d2
756  real(dp):: cijk
757  logical :: has_elastic3rd,has_elastic4th,has_elastic_displ
758 ! array
759 ! *************************************************************************
760 
761 !Reset output and flags
762  energy = zero
763  fcart = zero
764  strten = zero
765  has_elastic3rd    = .FALSE.
766  has_elastic4th    = .FALSE.
767  has_elastic_displ = .FALSE.
768  d1=0;d2=0
769 
770 !Set the flags
771  if(present(elastic3rd)) has_elastic3rd = .TRUE.
772  if(present(elastic4th)) then
773    has_elastic4th = .TRUE.
774    d1=1;d2=6
775  end if
776  if(present(elastic_displacement)) has_elastic_displ = .TRUE.
777 
778 !1-Treat 3rd order elastic constants
779  if (has_elastic3rd.or.has_elastic4th) then
780    do alpha=1,6
781      do beta=1,6
782        do gamma=1,6
783          cijk = ncell*elastic3rd(alpha,beta,gamma)
784 !        Accumulate energy
785          energy = energy + sixth*cijk*strain(alpha)*strain(beta)*strain(gamma)
786 !        Accumulate stresses contributions
787          strten(alpha)=strten(alpha)+ half*cijk*strain(beta)*strain(gamma)
788          do delta=d1,d2
789            cijk = ncell*elastic4th(alpha,beta,gamma,delta)
790 !          Accumulate energy
791            energy = energy + (1/24.)*cijk*strain(alpha)*strain(beta)*&
792 &                                                  strain(gamma)*strain(delta)
793 !          Accumulate stresses contributions
794            strten(alpha)=strten(alpha)+ sixth*cijk*strain(beta)*strain(gamma)*&
795 &                                                           strain(delta)
796          end do
797        end do
798      end do
799    end do
800  end if
801 
802 !2-Part due to the internat strain
803  if(has_elastic_displ)then
804    ii = 1
805    do ia = 1,natom
806      do mu = 1,3
807        do beta=1,6
808          do alpha=1,6
809            cijk = elastic_displacement(alpha,beta,mu,ii)
810 !          Accumulte for this atom
811            energy = energy + sixth*cijk*strain(alpha)*strain(beta)*disp(mu,ia)
812            fcart(mu,ia) = fcart(mu,ia)   +  half*cijk*strain(alpha)*strain(beta)
813            strten(alpha) = strten(alpha) +  half*cijk*strain(beta)*disp(mu,ia)
814          end do
815        end do
816      end do
817      ii = ii +1
818 !    Reset to 1 if the number of atoms is superior than in the initial cell
819      if(ii==natom_uc+1) ii = 1
820    end do
821  end if
822 
823 end subroutine anharmonics_terms_evaluateElastic