TABLE OF CONTENTS


ABINIT/m_harmonics_terms [ Functions ]

[ Top ] [ Functions ]

NAME

 m_harmonics_term

FUNCTION

 Module with datatype and tools for the harmonics 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_harmonics_terms
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_supercell,only: getPBCIndexes_supercell
32  use m_xmpi,only : xmpi_sum
33  use m_ifc
34 
35  implicit none
36 
37  public :: harmonics_terms_init
38  public :: harmonics_terms_free
39  public :: harmonics_terms_applySumRule
40  public :: harmonics_terms_evaluateIFC
41  public :: harmonics_terms_evaluateElastic
42  public :: harmonics_terms_setEffectiveCharges
43  public :: harmonics_terms_setDynmat
44  public :: harmonics_terms_setInternalStrain

defs_abitypes/harmonics_terms_type [ Types ]

[ Top ] [ defs_abitypes ] [ Types ]

NAME

 harmonics_terms_type

FUNCTION

 datatype for harmonic part of effective potential.

SOURCE

56  type, public :: harmonics_terms_type
57 
58    integer :: nqpt
59 !   Number of qpoints
60 
61    real(dp) :: epsilon_inf(3,3)
62 !   epsilon_inf(3,3)
63 !   Dielectric tensor
64 
65    real(dp) :: elastic_constants(6,6)
66 !   elastic_constant(6,6)
67 !   Elastic tensor Hartree
68 
69    real(dp), allocatable :: strain_coupling(:,:,:)
70 !   strain_coupling(6,3,natom)
71 !   internal strain tensor
72 
73    real(dp), allocatable :: zeff(:,:,:)
74 !   zeff(3,3,natom) Effective charges
75 
76    type(ifc_type) :: ifcs
77 !   type with ifcs constants (short + ewald)
78 !   also contains the number of cell and the indexes
79 
80    real(dp), allocatable :: qpoints(:,:)
81 !   qph1l(3,nqpt)
82 !   List of qpoints wavevectors
83 
84    real(dp), allocatable :: dynmat(:,:,:,:,:,:)
85 !   dynmat(2,3,natom,3,natom,nqpt)
86 !   dynamical matrix for each q points
87 
88    real(dp), allocatable :: phfrq(:,:)
89 !   phfrq(3*natom,nqpt)
90 !   array with all phonons frequencies for each q points in Hartree/cm
91 
92  end type harmonics_terms_type

m_harmonics_terms/harmonics_terms_evaluateElastic [ Functions ]

[ Top ] [ m_harmonics_terms ] [ Functions ]

NAME

  harmonics_terms_evaluateElastic

FUNCTION

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

INPUTS

  elastic_constants(6,6) = elastic constants in Hartree
  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
  ncell = total number of cell
  strain_coupling(6,3,natom) = internal strain coupling parameters
  strain(6) = strain between configuration and the reference

OUTPUT

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

PARENTS

      m_effective_potential

CHILDREN

      wrtout

SOURCE

724 subroutine harmonics_terms_evaluateElastic(elastic_constants,disp,energy,fcart,natom,natom_uc,ncell,&
725 &                                          strain_coupling,strten,strain)
726 
727 
728 !This section has been created automatically by the script Abilint (TD).
729 !Do not modify the following lines by hand.
730 #undef ABI_FUNC
731 #define ABI_FUNC 'harmonics_terms_evaluateElastic'
732 !End of the abilint section
733 
734  real(dp),intent(out):: energy
735  integer, intent(in) :: natom,natom_uc,ncell
736 ! array
737  real(dp),intent(in) :: elastic_constants(6,6),strain_coupling(6,3,natom)
738  real(dp),intent(out):: strten(6)
739  real(dp),intent(out):: fcart(3,natom)
740  real(dp),intent(in) :: disp(3,natom)
741  real(dp),intent(in) :: strain(6)
742 
743  !Local variables-------------------------------
744 ! scalar
745  integer :: ia,ii,mu,alpha,beta
746  real(dp):: cij
747 ! array
748 ! *************************************************************************
749 
750  energy = zero
751  fcart = zero
752  strten = zero
753 
754 !1- Part due to elastic constants
755  do alpha=1,6
756    do beta=1,6
757      cij = ncell*elastic_constants(alpha,beta)
758      energy = energy + half*cij*strain(alpha)*strain(beta)
759      strten(alpha) = strten(alpha) + cij*strain(beta)
760    end do
761  end do
762 
763 !2-Part due to the internal strain coupling parameters
764  ii = 1
765  do ia = 1,natom
766    do mu = 1,3
767      do alpha=1,6
768        cij = strain_coupling(alpha,mu,ii)
769 !      Accumulte for this atom
770        energy = energy + half*cij*strain(alpha)*disp(mu,ia)
771        fcart(mu,ia)  = fcart(mu,ia)  + half*cij*strain(alpha)
772        strten(alpha) = strten(alpha) + half*cij*disp(mu,ia)
773      end do
774    end do
775    ii = ii +1
776 !  Reset to 1 if the number of atoms is superior than in the initial cell
777    if(ii==natom_uc+1) ii = 1
778  end do
779 
780 end subroutine  harmonics_terms_evaluateElastic

m_harmonics_terms/harmonics_terms_evaluateIFC [ Functions ]

[ Top ] [ m_harmonics_terms ] [ Functions ]

NAME

  harmonics_terms_evaluateIFC

FUNCTION

  This fonction compute the contribution of the ifc harmonic part of
  the energy and forces.

INPUTS

  atmfrc(3,natom_uc,3,natom_uc,nrpt) = atomic force constants
  disp(3,natom_sc) = atomics displacement between configuration and the reference
  ncell = total number of cell to treat
  nrpt  = total number of rpt to treat
  natom_sc = number of atoms in the supercell
  natom_uc = number of atoms in the unit cell
  nrpt  = number of rpt
  atmrpt_index(nrpt,cell) = For each cell in the supercell and each rpt,
                            give the index of the first atoms in the rpt cell
  rpt(nrpt) = index of rpt in  atmfrc (6th dimension)
  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

 PARENT
   effective_potential_evaluate

PARENTS

      m_effective_potential

CHILDREN

      wrtout

SOURCE

608 subroutine harmonics_terms_evaluateIFC(atmfrc,disp,energy,fcart,natom_sc,natom_uc,ncell,nrpt,&
609 &                                      atmrpt_index,index_cells,sc_size,rpt,comm)
610 
611 
612 !This section has been created automatically by the script Abilint (TD).
613 !Do not modify the following lines by hand.
614 #undef ABI_FUNC
615 #define ABI_FUNC 'harmonics_terms_evaluateIFC'
616 !End of the abilint section
617 
618  implicit none
619 
620 !Arguments -------------------------------
621 ! scalars
622   real(dp),intent(out) :: energy
623   integer,intent(in) :: natom_uc,natom_sc,ncell,nrpt
624   integer,intent(in) :: comm
625 ! array
626   integer,intent(in) :: sc_size(3),atmrpt_index(nrpt,ncell)
627   integer,intent(in) ::  index_cells(4,ncell),rpt(nrpt)
628   real(dp),intent(in) :: atmfrc(3,natom_uc,3,natom_uc,nrpt)
629   real(dp),intent(in) :: disp(3,natom_sc)
630   real(dp),intent(out) :: fcart(3,natom_sc)
631 
632 !Local variables-------------------------------
633 ! scalar
634   integer :: i1,i2,i3,ia,ib,icell,ierr,irpt,irpt_tmp,ii,jj,kk,ll
635   integer :: mu,nu
636   real(dp):: disp1,disp2,ifc,tmp,tmp2
637 ! array
638   character(500) :: msg
639 
640 ! *************************************************************************
641 
642   if (any(sc_size <= 0)) then
643     write(msg,'(a,a)')' sc_size can not be inferior or equal to zero'
644     MSG_ERROR(msg)
645   end if
646 
647 ! Initialisation of variables
648   energy   = zero
649   fcart(:,:) = zero
650 
651   do icell = 1,ncell
652     i1 = index_cells(1,icell)
653     i2 = index_cells(2,icell)
654     i3 = index_cells(3,icell)
655 !   index of the first atom in the current cell
656     ii = index_cells(4,icell)
657     do irpt_tmp = 1,nrpt
658       irpt = rpt(irpt_tmp)
659 !     index of the first atom in the irpt cell
660       jj = atmrpt_index(irpt_tmp,icell)
661 !     Loop over the atom in the cell
662       do ib = 1, natom_uc
663         ll = jj + ib
664         do nu=1,3
665           disp2 = disp(nu,ll)
666           do ia = 1, natom_uc
667             kk = ii + ia
668             do mu=1,3
669               disp1 = disp(mu,kk)
670               ifc = atmfrc(mu,ia,nu,ib,irpt)
671 !              if(abs(ifc) > tol10)then
672                 tmp = disp2 * ifc
673 !               accumule energy
674                 tmp2 = disp1*tmp
675                 energy =  energy + tmp2
676 !               accumule forces
677                 fcart(mu,kk) = fcart(mu,kk) + tmp
678 !              end if
679             end do
680           end do
681         end do
682       end do
683     end do
684   end do
685 
686   energy = half * energy
687 
688 ! MPI_SUM
689   call xmpi_sum(energy, comm, ierr)
690   call xmpi_sum(fcart , comm, ierr)
691 
692 end subroutine harmonics_terms_evaluateIFC

m_harmonics_terms/harmonics_terms_init [ Functions ]

[ Top ] [ m_harmonics_terms ] [ Functions ]

NAME

 harmonics_terms_init

FUNCTION

 Initialize harmonics_terms datatype

INPUTS

 ifc<type(ifc_type)> = interatomic forces constants
 natom = number of atoms in primitive cell
 nrpt = number rpt (cell) in the ifc
 dynmat(2,3,natom,3,natom,3,nqpt) = optional, dynamical matricies for each q-point
 epsilon_inf(3,3) = optional, dielectric tensor
 elastic_constant(6,6) = optional, elastic constant
 strain_coupling(6,3,natom) = optional, internal strain coupling parameters
 nqpt = optional, number of q-points
 phfrq(3*natom,nqpt) = optional,phonons frequencies for each q points in Hartree/cm
 qpoints(3,nqpt) = list of qpoints wavevectors
 zeff(3,3,natom) = optional,effective charges

OUTPUT

 harmonics_terms<type(harmonics_terms_type)> = harmonics_terms datatype to be initialized

PARENTS

      m_effective_potential

CHILDREN

      wrtout

SOURCE

130 subroutine harmonics_terms_init(harmonics_terms,ifcs,natom,nrpt,&
131 &                               dynmat,epsilon_inf,elastic_constants,strain_coupling,&
132 &                               nqpt,phfrq,qpoints,zeff)
133 
134 
135 !This section has been created automatically by the script Abilint (TD).
136 !Do not modify the following lines by hand.
137 #undef ABI_FUNC
138 #define ABI_FUNC 'harmonics_terms_init'
139 !End of the abilint section
140 
141  implicit none
142 
143 !Arguments ------------------------------------
144 !scalars
145  integer, intent(in) :: natom,nrpt
146 !arrays
147  type(ifc_type),intent(in) :: ifcs
148  type(harmonics_terms_type), intent(out) :: harmonics_terms
149  integer, optional,intent(in) :: nqpt
150  real(dp),optional,intent(in) :: epsilon_inf(3,3),dynmat(:,:,:,:,:,:)
151  real(dp),optional,intent(in) :: elastic_constants(6,6)
152  real(dp),optional,intent(in) :: strain_coupling(6,3,natom),zeff(3,3,natom)
153  real(dp),optional,intent(in) :: phfrq(:,:),qpoints(:,:)
154 !Local variables-------------------------------
155 !scalar
156 !arrays
157  character(len=500) :: msg
158 
159 ! *************************************************************************
160 
161  call harmonics_terms_free(harmonics_terms)
162 
163 ! Do some Checks
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  if (nrpt < 1) then
172    write(msg, '(a,a,a,i10,a)' )&
173 &   'The cell must have at least one rpt point.',ch10,&
174 &   'The number of rpt points is  ',nrpt,'.'
175    MSG_BUG(msg)
176  end if
177 
178  if (nrpt /= ifcs%nrpt) then
179    write(msg, '(3a,i5,a,i5,a)' )&
180 &   'nrpt must have the same dimension as ifcs.',ch10,&
181 &   'The number of cell is  ',nrpt,' instead of ',ifcs%nrpt,'.'
182    MSG_BUG(msg)
183  end if
184 
185  if(present(nqpt).and.(.not.present(dynmat).or.&
186 &                      .not.present(qpoints)   .or.&
187 &                      .not.present(phfrq)))then
188    write(msg, '(a)' )&
189 &   'nqpt is specified but dynamt,qpoints or phfrq are not.'
190    MSG_BUG(msg)
191  end if
192 
193  if(.not.present(nqpt).and.(present(dynmat).or.&
194 &                      present(qpoints)   .or.&
195 &                      present(phfrq)))then
196    write(msg, '(a)' )&
197 &   ' dynamt,qpoints or phfrq are specified but nqpt is not.'
198    MSG_BUG(msg)
199  end if
200 
201 !Set number of cell
202  harmonics_terms%ifcs%nrpt = nrpt
203 
204 !Allocation of total ifc
205  ABI_ALLOCATE(harmonics_terms%ifcs%atmfrc,(3,natom,3,natom,nrpt))
206  harmonics_terms%ifcs%atmfrc(:,:,:,:,:) = ifcs%atmfrc(:,:,:,:,:)
207 
208 !Allocation of ewald part of ifc
209  ABI_ALLOCATE(harmonics_terms%ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt))
210  harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,:) = ifcs%ewald_atmfrc(:,:,:,:,:)
211 
212 !Allocation of short range part of ifc
213  ABI_ALLOCATE(harmonics_terms%ifcs%short_atmfrc,(3,natom,3,natom,nrpt))
214  harmonics_terms%ifcs%short_atmfrc(:,:,:,:,:) = ifcs%short_atmfrc(:,:,:,:,:)
215 
216 !Allocation of cell of ifc
217  ABI_ALLOCATE(harmonics_terms%ifcs%cell,(3,nrpt))
218  harmonics_terms%ifcs%cell(:,:) = ifcs%cell(:,:)
219 
220 !Allocation of the dynamical matrix
221  harmonics_terms%nqpt = 0
222  if(present(nqpt).and.present(dynmat).and.present(qpoints).and.present(phfrq))then
223    call harmonics_terms_setDynmat(dynmat,harmonics_terms,natom,nqpt,&
224 &                                 harmonics_terms%phfrq,harmonics_terms%qpoints)
225  end if
226 
227 !Allocation of the elastic constants
228  harmonics_terms%elastic_constants = zero
229  if (present(elastic_constants)) then
230    harmonics_terms%elastic_constants = elastic_constants
231  end if
232 
233 !Allication of the dielectric tensor
234  harmonics_terms%epsilon_inf = zero
235  if (present(epsilon_inf)) then
236    harmonics_terms%epsilon_inf = epsilon_inf
237  end if
238 
239 !Allocation of Effective charges array
240  ABI_ALLOCATE(harmonics_terms%zeff,(3,3,natom))
241  harmonics_terms%zeff = zero
242  if (present(zeff)) then
243    call harmonics_terms_setEffectiveCharges(harmonics_terms,natom,zeff)
244    harmonics_terms%zeff = zeff
245  end if
246 
247 !Allocation of internal strain tensor
248  ABI_ALLOCATE(harmonics_terms%strain_coupling,(6,3,natom))
249  harmonics_terms%strain_coupling = zero
250  if (present(strain_coupling)) then
251    call harmonics_terms_setInternalStrain(harmonics_terms,natom,strain_coupling)
252  end if
253 
254 end subroutine harmonics_terms_init