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

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

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

SOURCE

723 subroutine anharmonics_terms_evaluateIFCStrainCoupling(phonon_strain,disp,energy,fcart,natom,natom_uc,&
724 &                                                      sc_size,strain,strten,cells,ncell,&
725 &                                                      index_cells,comm)
726 
727 !Arguments -------------------------------
728 ! scalars
729   real(dp),intent(out) :: energy
730   integer,intent(in) :: natom,natom_uc,ncell
731   integer,intent(in) :: comm
732 ! array
733   integer,intent(in) ::   cells(ncell),index_cells(ncell,3)
734   integer,intent(in) :: sc_size(3)
735   type(ifc_type),intent(in) :: phonon_strain(6)
736   real(dp),intent(in) :: disp(3,natom)
737   real(dp),intent(out) :: fcart(3,natom)
738   real(dp),intent(out) :: strten(6)
739   real(dp),intent(in) :: strain(6)
740 !Local variables-------------------------------
741 ! scalar
742   integer :: alpha
743   integer :: i1,i2,i3,ia,ib,icell,ii
744   integer :: irpt,jj,kk,ll,mu,nu
745   integer :: ierr
746   real(dp):: ifc
747 ! array
748   integer :: cell_atom2(3)
749   character(500) :: msg
750 
751 ! *************************************************************************
752 
753   if (any(sc_size <= 0)) then
754     write(msg,'(a,a)')' sc_size can not be inferior or equal to zero'
755     ABI_ERROR(msg)
756   end if
757 
758 ! Initialisation of variables
759   energy   = zero
760   fcart(:,:) = zero
761   strten(:) = zero
762 
763   do icell = 1,ncell
764     ii = (cells(icell)-1)*natom_uc
765     i1=index_cells(icell,1); i2=index_cells(icell,2); i3=index_cells(icell,3)
766     do alpha=1,6
767       do irpt = 1,phonon_strain(alpha)%nrpt
768 !       get the cell of atom2  (0 0 0, 0 0 1...)
769         cell_atom2(1) =  i1 + phonon_strain(alpha)%cell(1,irpt)
770         cell_atom2(2) =  i2 + phonon_strain(alpha)%cell(2,irpt)
771         cell_atom2(3) =  i3 + phonon_strain(alpha)%cell(3,irpt)
772         call getPBCIndexes_supercell(cell_atom2(1:3),sc_size(1:3))
773 !       index of the second atom in the displacement array
774         jj = (cell_atom2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
775 &            (cell_atom2(2)-1)*sc_size(3)*natom_uc+&
776 &            (cell_atom2(3)-1)*natom_uc
777         do ib = 1, natom_uc
778           ll = jj + ib
779           do nu=1,3
780             do ia = 1, natom_uc
781               kk = ii + ia
782               do mu=1,3
783                 ifc = phonon_strain(alpha)%atmfrc(mu,ia,nu,ib,irpt)
784 !               accumule energy
785                 energy =  energy + sixth*strain(alpha)*disp(mu,kk)*disp(nu,ll)*ifc
786 !               accumule forces
787                 fcart(mu,kk) = fcart(mu,kk) + half*strain(alpha)*disp(nu,ll)*ifc
788 !               accumule stresses
789                 strten(alpha) = strten(alpha) + half*disp(mu,kk)*disp(nu,ll)*ifc
790               end do
791             end do
792           end do
793         end do
794       end do
795     end do
796   end do
797 
798 ! MPI_SUM
799   call xmpi_sum(energy, comm, ierr)
800   call xmpi_sum(fcart , comm, ierr)
801   call xmpi_sum(strten, comm, ierr)
802 
803 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

SOURCE

124 subroutine anharmonics_terms_init(anharmonics_terms,natom,ncoeff,&
125 &                                 bounded,elastic3rd,elastic4th,elastic_displacement,&
126 &                                 phonon_strain,coeffs)
127 
128 !Arguments ------------------------------------
129 !scalars
130  integer, intent(in) :: natom,ncoeff
131  type(anharmonics_terms_type), intent(out) :: anharmonics_terms
132  real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
133  real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
134  type(polynomial_coeff_type),optional :: coeffs(ncoeff)
135  type(ifc_type),optional,intent(in) :: phonon_strain(6)
136  logical,optional,intent(in) :: bounded
137 !arrays
138 !Local variables-------------------------------
139 !scalar
140 !arrays
141  character(len=500) :: msg
142 
143 ! *************************************************************************
144 
145  call anharmonics_terms_free(anharmonics_terms)
146 
147 ! Check the number of atoms
148  if (natom < 1) then
149    write(msg, '(a,a,a,i10,a)' )&
150 &   'The cell must have at least one atom.',ch10,&
151 &   'The number of atom is  ',natom,'.'
152    ABI_BUG(msg)
153  end if
154 
155 !Allocation of phonon strain coupling array (3rd order)
156  if(present(phonon_strain)) then
157    call anharmonics_terms_setStrainPhononCoupling(anharmonics_terms,natom,phonon_strain)
158  end if
159 
160 !Set the 3rd order elastic tensor
161  anharmonics_terms%elastic3rd = zero
162  anharmonics_terms%has_elastic3rd = .FALSE.
163  if(present(elastic3rd))then
164    call anharmonics_terms_setElastic3rd(anharmonics_terms,elastic3rd)
165  end if
166 
167 !Set the 3rd order elastic tensor
168  anharmonics_terms%elastic4th = zero
169  anharmonics_terms%has_elastic4th = .FALSE.
170  if(present(elastic4th))then
171    call anharmonics_terms_setElastic4th(anharmonics_terms,elastic4th)
172  end if
173 
174 !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement
175  if(present(elastic_displacement))then
176    call anharmonics_terms_setElasticDispCoupling(anharmonics_terms,natom,elastic_displacement)
177  end if
178 
179  anharmonics_terms%ncoeff = 0
180 
181 !Allocation of the coefficient
182   if(present(coeffs))then
183    if(ncoeff /= size(coeffs))then
184      write(msg, '(a)' )&
185 &        ' ncoeff has not the same size than coeffs array, '
186      ABI_BUG(msg)
187    end if
188    call anharmonics_terms_setCoeffs(coeffs,anharmonics_terms,ncoeff)
189  end if
190 
191 !Set the flag bounded
192  if(present(bounded))then
193    anharmonics_terms%bounded = bounded
194  else
195    anharmonics_terms%bounded = .FALSE.
196  end if
197 
198 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

SOURCE

609 subroutine anharmonics_terms_evaluateElastic(disp,energy,fcart,natom,natom_uc,ncell,strten,strain,&
610 &                                            elastic3rd,elastic4th,elastic_displacement)
611 
612  real(dp),intent(out):: energy
613  integer, intent(in) :: natom,natom_uc,ncell
614 ! array
615  real(dp),optional,intent(in) :: elastic3rd(6,6,6),elastic4th(6,6,6,6)
616  real(dp),optional,intent(in) :: elastic_displacement(6,6,3,natom)
617  real(dp),intent(out):: strten(6)
618  real(dp),intent(out):: fcart(3,natom)
619  real(dp),intent(in) :: disp(3,natom)
620  real(dp),intent(in) :: strain(6)
621 
622  !Local variables-------------------------------
623 ! scalar
624  integer :: ia,ii,mu,alpha,beta,gamma,delta,d1,d2
625  real(dp):: cijk
626  logical :: has_elastic3rd,has_elastic4th,has_elastic_displ
627 ! array
628 ! *************************************************************************
629 
630 !Reset output and flags
631  energy = zero
632  fcart = zero
633  strten = zero
634  has_elastic3rd    = .FALSE.
635  has_elastic4th    = .FALSE.
636  has_elastic_displ = .FALSE.
637  d1=0;d2=0
638 
639 !Set the flags
640  if(present(elastic3rd)) has_elastic3rd = .TRUE.
641  if(present(elastic4th)) then
642    has_elastic4th = .TRUE.
643    d1=1;d2=6
644  end if
645  if(present(elastic_displacement)) has_elastic_displ = .TRUE.
646 
647 !1-Treat 3rd order elastic constants
648  if (has_elastic3rd.or.has_elastic4th) then
649    do alpha=1,6
650      do beta=1,6
651        do gamma=1,6
652          cijk = ncell*elastic3rd(alpha,beta,gamma)
653 !        Accumulate energy
654          energy = energy + sixth*cijk*strain(alpha)*strain(beta)*strain(gamma)
655 !        Accumulate stresses contributions
656          strten(alpha)=strten(alpha)+ half*cijk*strain(beta)*strain(gamma)
657          do delta=d1,d2
658            cijk = ncell*elastic4th(alpha,beta,gamma,delta)
659 !          Accumulate energy
660            energy = energy + (1/24.)*cijk*strain(alpha)*strain(beta)*&
661 &                                                  strain(gamma)*strain(delta)
662 !          Accumulate stresses contributions
663            strten(alpha)=strten(alpha)+ sixth*cijk*strain(beta)*strain(gamma)*&
664 &                                                           strain(delta)
665          end do
666        end do
667      end do
668    end do
669  end if
670 
671 !2-Part due to the internat strain
672  if(has_elastic_displ)then
673    ii = 1
674    do ia = 1,natom
675      do mu = 1,3
676        do beta=1,6
677          do alpha=1,6
678            cijk = elastic_displacement(alpha,beta,mu,ii)
679 !          Accumulte for this atom
680            energy = energy + sixth*cijk*strain(alpha)*strain(beta)*disp(mu,ia)
681            fcart(mu,ia) = fcart(mu,ia)   +  half*cijk*strain(alpha)*strain(beta)
682            strten(alpha) = strten(alpha) +  half*cijk*strain(beta)*disp(mu,ia)
683          end do
684        end do
685      end do
686      ii = ii +1
687 !    Reset to 1 if the number of atoms is superior than in the initial cell
688      if(ii==natom_uc+1) ii = 1
689    end do
690  end if
691 
692 end subroutine anharmonics_terms_evaluateElastic