TABLE OF CONTENTS


m_polynomial_coeff/check_irreducibility [ Functions ]

[ Top ] [ Functions ]

NAME

 check_irreducibility

FUNCTION

 checks irreducibility of combination of terms with respect to crystal symmetry

INPUTS

 combination: combination of integers icoeff representing term A_x-B_x to check
 list_combination: list of combination against which combinatino will be checked
 list_symcoeff: list of coefficients containig information connecting the integers in combination
                and list_combination to bodies (A_x-B_x) of difference of atomic displacements
 list_symstr: list of strain with symmetry properties connecting icoeff>ncoeff_sym to strains of given direction
 ncoeff_sym:  number of bodies (A_x-B_x) containing all symmetric equivalents
 nsym:        number of symmetries in the crystal system
 ncombination: number of combinations -> defines size of list_combination
 ndisp:        number of displacements that is maximum order of combinations
 index_irred:  list of index of irreducible terms (that are non-zero) in list_combination

OUTPUT

 logical :: irreducible -> TRUE: no other equal term exists in list_combination
                        -> FALSE: a other equivalent term exists already

 SIDE_EFFECT: if combination is irreducible (that means irreducible = .FALSE.) the index of the irreduc!!              irreducible coefficient is added to index_irred

SOURCE

4363 function check_irreducibility(combination,list_combination,list_symcoeff,list_symstr,ncoeff_sym,nsym,&
4364 &                             ncombination,ndisp,index_irred)  result(irreducible)
4365 !Arguments ------------------------------------
4366  implicit none
4367 
4368 !Arguments ------------------------------------
4369  !scalar
4370  integer,intent(in) :: ncoeff_sym,ncombination,ndisp,nsym
4371  integer,intent(inout),allocatable :: index_irred(:)
4372  logical :: irreducible !output
4373  !array
4374  integer,intent(inout) :: combination(ndisp),list_combination(ndisp,ncombination)
4375  integer,intent(in) :: list_symcoeff(6,ncoeff_sym,nsym)
4376  integer,intent(in) :: list_symstr(6,nsym,2)
4377 !local
4378 !variable
4379   integer :: icombi, isym
4380 !array
4381   integer :: combination_cmp_tmp(ndisp)
4382   integer,allocatable :: index_irred_tmp(:)
4383 ! *************************************************************************
4384  ! sort input combinations
4385  !write(std_out,*) "DEBUG call sort_combination_lsit"
4386  call sort_combination_list(list_combination,size(list_combination,1),size(list_combination,2))
4387  !write(std_out,*) "DEBUG sort_combination in check_irreducibility"
4388  call sort_combination(combination,size(combination))
4389 
4390  !Initialize variables
4391  icombi = 1
4392  irreducible = .TRUE.
4393  !do icombi=1,ncombination
4394  ! write(std_out,*) "DEBUG: list_combination i:",icombi,": ", list_combination(:,icombi)
4395  !enddo
4396  !MS DEBUG
4397  ! write(std_out,*) "DEBUG: shape(list_combination)", shape(list_combination(:,:))
4398  !write(std_out,*) "DEBUG: combination: ", combination
4399  !write(std_out,*) "DEBUG: ncombination: ", ncombination
4400 
4401  !Loop over non reducable combinations
4402  do while(icombi <= size(index_irred) .and. irreducible .eqv. .TRUE.)
4403    !if(all(list_combination(:,index_irred(icombi)) /= 0))then !.and. any(list_combination(:,i) <= ncoeff_symsym))then
4404      !If term j is equivalent to term i delete it
4405      if(all(list_combination(:,index_irred(icombi)) == combination(:)))then
4406         irreducible = .FALSE.
4407         return
4408      else !else loop over symmetries to find symmetry operation that projects term j on i
4409        isym = 2
4410        do while(isym <= nsym)
4411           !Get equivalent term indexes for symmetry isym
4412          call symcomb(combination, combination_cmp_tmp, isym)
4413          !Sort the new symmetric indexes for comparision
4414          call sort_combination(combination_cmp_tmp,size(combination_cmp_tmp))
4415          !Compare. If equivalent break (term not irreducible
4416           if(all(list_combination(:,index_irred(icombi)) == combination_cmp_tmp(:)))then
4417              irreducible = .FALSE.
4418              return
4419           else
4420              isym = isym + 1
4421           endif
4422        enddo !isym
4423      endif ! all(list_combinaton...)
4424    !end if ! any(list_combination /= 0 )
4425    icombi = icombi + 1
4426 end do !i=1,ncombination
4427 
4428 !If no equivalent found put ncombination +1 to index of irreducible combis
4429 ABI_MALLOC(index_irred_tmp,(size(index_irred)))
4430 index_irred_tmp = index_irred
4431 ABI_FREE(index_irred)
4432 ABI_MALLOC(index_irred,(size(index_irred_tmp)+1))
4433 index_irred(:size(index_irred_tmp)) = index_irred_tmp
4434 ABI_FREE(index_irred_tmp)
4435 index_irred(size(index_irred)) = ncombination + 1
4436 !write(std_out,*) "DEBUG: index_irred", index_irred
4437 return
4438 
4439 contains
4440   subroutine symcomb(combination, combination_cmp_tmp, isym)
4441     integer, intent(inout) :: combination(:), combination_cmp_tmp(:), isym
4442     integer :: idisp, istrain
4443     do idisp=1,ndisp
4444       if(combination(idisp) /= 0 .and. combination(idisp) <= ncoeff_sym)then
4445         combination_cmp_tmp(idisp)=list_symcoeff(6,combination(idisp),isym)
4446       else if(combination(idisp) > ncoeff_sym)then
4447         istrain = combination(idisp) - ncoeff_sym
4448         combination_cmp_tmp(idisp)=list_symstr(istrain,isym,1) + ncoeff_sym
4449       else
4450         combination_cmp_tmp(idisp) = 0
4451       endif
4452     enddo
4453   end subroutine symcomb
4454 
4455 end function check_irreducibility

m_polynomial_coeff/coeffs_compare [ Functions ]

[ Top ] [ Functions ]

NAME

  equal

FUNCTION

INPUTS

OUTPUT

SOURCE

4026 function coeffs_compare(c1,c2) result (res)
4027 !Arguments ------------------------------------
4028  implicit none
4029 
4030 !Arguments ------------------------------------
4031   type(polynomial_coeff_type), intent(in) :: c1,c2
4032   logical :: res
4033 !local
4034 !variable
4035   integer :: iterm1,iterm2
4036 !array
4037   !integer,allocatable :: blkval(:,:)
4038   integer :: blkval(2, max(c1%nterm,c2%nterm))
4039 ! *************************************************************************
4040   res = .false.
4041   blkval = 0
4042   do iterm1=1,c1%nterm
4043     if(blkval(1,iterm1)==1)cycle!already found
4044     do iterm2=1,c2%nterm
4045       if(blkval(2,iterm2)==1)cycle!already found
4046       if(c1%terms(iterm1)==c2%terms(iterm2)) then
4047         blkval(1,iterm1) = 1
4048         blkval(2,iterm2) = 1
4049       end if
4050     end do
4051   end do
4052   if(.not.any(blkval(:,:)==0))res = .true.
4053 
4054 end function coeffs_compare

m_polynomial_coeff/coeffs_list_append [ Functions ]

[ Top ] [ Functions ]

NAME

 coeff_list_append

FUNCTION

 append on coeff to the list of coeffs: coeff_list

INPUTS

OUTPUT

SOURCE

4167 subroutine coeffs_list_append(coeff_list,coeff, check)
4168 !Arguments ------------------------------------
4169  implicit none
4170 
4171 !Arguments ------------------------------------
4172  type(polynomial_coeff_type), allocatable, intent(inout) ::coeff_list(:)
4173  type(polynomial_coeff_type), intent(inout) :: coeff
4174  logical, intent(in) :: check
4175  !local variable
4176  type(polynomial_coeff_type) :: tmp(size(coeff_list))
4177  integer :: n1, n2
4178  !array
4179  ! *************************************************************************
4180 
4181  n1=size(coeff_list)
4182  n2= n1+1
4183 
4184  ! copy list1 to tmp
4185  tmp=coeff_list
4186 
4187  ! allocate new list1
4188  call polynomial_coeff_list_free(coeff_list)
4189  ABI_MALLOC(coeff_list,(n2))
4190 
4191  call coeffs_list_copy(coeff_list, tmp)
4192  call polynomial_coeff_init(coeff%coefficient,coeff%nterm,coeff_list(n2),coeff%terms,&
4193    &                                 coeff%name,check=check)
4194 end subroutine coeffs_list_append

m_polynomial_coeff/coeffs_list_conc [ Functions ]

[ Top ] [ Functions ]

NAME

 coeff_list_conc

FUNCTION

 Concatenate list1 and list2 of type polynomial_coeff and store it in list_out

INPUTS

OUTPUT

SOURCE

4072 function coeffs_list_conc(coeff_list1,coeff_list2) result (coeff_list_out)
4073 !Arguments ------------------------------------
4074  implicit none
4075 
4076 !Arguments ------------------------------------
4077   type(polynomial_coeff_type), intent(in) :: coeff_list1(:),coeff_list2(:)
4078   type(polynomial_coeff_type) :: coeff_list_out(size(coeff_list1)+size(coeff_list2))
4079 !local
4080 !variable
4081   integer :: ncoeff1,ncoeff2,ncoeff_out,i,j
4082 !array
4083 ! *************************************************************************
4084 
4085 !Get sizes of coeff_list1/2
4086  ncoeff1 = size(coeff_list1)
4087  ncoeff2 = size(coeff_list2)
4088  ncoeff_out = ncoeff1 + ncoeff2
4089 
4090 ! ABI_MALLOC(coeff_list_out,(ncoeff_out))
4091  do i=1,ncoeff_out
4092     if(i<=ncoeff1)then
4093        call polynomial_coeff_init(coeff_list1(i)%coefficient,coeff_list1(i)%nterm,coeff_list_out(i),coeff_list1(i)%terms,&
4094 &                                 coeff_list1(i)%name,check=.TRUE.)
4095     else
4096        j=i-ncoeff1
4097        call polynomial_coeff_init(coeff_list2(j)%coefficient,coeff_list2(j)%nterm,coeff_list_out(i),coeff_list2(j)%terms,&
4098 &                                 coeff_list2(j)%name,check=.TRUE.)
4099     endif
4100  enddo
4101 
4102 end function coeffs_list_conc

m_polynomial_coeff/coeffs_list_conc_onsite [ Functions ]

[ Top ] [ Functions ]

NAME

 coeff_list_conc_onsite

FUNCTION

 Concatenate list1 and list2 of type polynomial_coeff and store it in list1

INPUTS

OUTPUT

SOURCE

4119 subroutine coeffs_list_conc_onsite(coeff_list1,coeff_list2)
4120 !Arguments ------------------------------------
4121  implicit none
4122 
4123 !Arguments ------------------------------------
4124  type(polynomial_coeff_type), allocatable, intent(inout) :: coeff_list1(:)
4125  type(polynomial_coeff_type), intent(in) ::coeff_list2(:)
4126  !local
4127  !variable
4128  integer :: ncoeff1,ncoeff2,ncoeff_out
4129  !array
4130  type(polynomial_coeff_type), allocatable :: coeff_list_tmp(:)
4131  ! *************************************************************************
4132 
4133  ncoeff1=size(coeff_list1)
4134  ncoeff2=size(coeff_list2)
4135  ncoeff_out = ncoeff1+ncoeff2
4136 
4137  ! copy list1 to tmp
4138  ABI_MALLOC(coeff_list_tmp,(ncoeff1))
4139  coeff_list_tmp=coeff_list1
4140 
4141  ! allocate new list1
4142  call polynomial_coeff_list_free(coeff_list1)
4143  ABI_MALLOC(coeff_list1,(ncoeff_out))
4144  !coeff_list1=coeff_list_tmp + coeff_list2
4145  !coeff_list1 = coeffs_list_conc(coeff_list_tmp, coeff_list2)
4146  coeff_list1(1:ncoeff1) = coeff_list_tmp
4147  coeff_list1(ncoeff1+1:ncoeff_out) = coeff_list2
4148  call polynomial_coeff_list_free(coeff_list_tmp)
4149 
4150 end subroutine coeffs_list_conc_onsite

m_polynomial_coeff/coeffs_list_copy [ Functions ]

[ Top ] [ Functions ]

NAME

 coeff_list_copy

FUNCTION

 Copy list1 to list2 of type polynomial_coeff

INPUTS

OUTPUT

SOURCE

4211 subroutine coeffs_list_copy(coeff_list_out,coeff_list_in)
4212 !Arguments ------------------------------------
4213  implicit none
4214 
4215 !Arguments ------------------------------------
4216   type(polynomial_coeff_type), intent(in) :: coeff_list_in(:)
4217   type(polynomial_coeff_type),intent(out) :: coeff_list_out(:)
4218 !local
4219 !variable
4220   integer :: ncoeff_in,ncoeff_out,ii
4221   logical :: check
4222   character(len=500):: message
4223 !array
4224 ! *************************************************************************
4225 
4226  check = .true.
4227 
4228 !Get size of coeff_lists
4229  ncoeff_in = size(coeff_list_in)
4230  ncoeff_out = size(coeff_list_out)
4231 
4232  if(ncoeff_in > ncoeff_out)then
4233    write(message,'(a,a,a)')'The input list of polynomial_coefficients is larger',ch10,&
4234 &                          'than the output list you want it assign to. Check size of lists.'
4235    ABI_ERROR(message)
4236  endif
4237 
4238  !Copy input list into output lists
4239  do ii=1,ncoeff_in
4240     call polynomial_coeff_init(coeff_list_in(ii)%coefficient,coeff_list_in(ii)%nterm,&
4241 &                             coeff_list_out(ii),coeff_list_in(ii)%terms,coeff_list_in(ii)%name,check)
4242  enddo
4243 end subroutine coeffs_list_copy

m_polynomial_coeff/computeCombinationFromList [ Functions ]

[ Top ] [ Functions ]

NAME

 computeCombinationFromList

FUNCTION

 Recursive routine to compute the order N of a all the possible coefficient
 from the list list_symcoeff and list_symstr.

INPUTS

 cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
 compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1
 list_coeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 list_str(nstr_sym,nsym) = array with the list of the strain  and the symmetrics
 index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0)
 icoeff = current indexes of the combination (start with 1)
 max_power_strain = maximum order of the strain of the strain phonon coupling
 nmodel_tot = current number of combination already computed (start with 0)
 natom = number of atoms in the unit cell
 ncoeff = number of coefficient for related to the atomic displacment into list_symcoeff
 nstr = number of coefficient for related to the strain into list_symstr
 nmodel = number of maximum models
 nrpt = number of cell
 nsym = number of symmetries in the system
     For example, the sum of all the term like (Sr_y-O_y)^odd, are 0 by symetrie in cubic system.
     Here, we build a list with: 0 this term is not allowed for odd
                                 1 this term is allowed for odd
 power_disp = initial power_disp to be computed (can be < power_disp_min,
              this routine will skip the first power_disp)
 power_disp_min = minimal power_disp to be computed
 power_disp_max = maximum power_disp to be computed
 symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
 nbody = optional, number of body for the coefficients, for example:
                   0 => all the terms
                   1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp  ...
 compute = logical, optional: TRUE if we store the coefficients
                              FALSE just to count the number of coefficient
 anharmstr = logical, optional : TRUE, the anharmonic strain are computed
                                   FALSE, (default) the anharmonic strain are not computed
 distributed = logical, optional : True, the coefficients will be distributed on the CPU
 only_odd_power = logical, optional : if TRUE return only odd power
 only_even_power= logical, optional : if TRUe return only even power

OUTPUT

 icoeff = current indexes of the cofficients (start we 1)
 nmodel_tot = current number of coefficients already computed (start we 0)
 list_combination = list of the possible combination of coefficients

SOURCE

3131 recursive subroutine computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
3132 &                                  index_coeff_in,list_combination,icoeff,max_power_strain,&
3133 &                                  natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp,power_disp_min,&
3134 &                                  power_disp_max,symbols,comm,nbody,only_odd_power,only_even_power,&
3135 &                                  compute,anharmstr,spcoupling,disp)
3136 
3137  implicit none
3138 
3139 !Arguments ---------------------------------------------
3140 !scalar
3141  integer,intent(in) :: natom,ncoeff,ncoeff_sym,power_disp,power_disp_min,power_disp_max
3142  integer,intent(in) :: max_power_strain,nmodel,nsym,nrpt,nstr,comm,icoeff
3143  integer,intent(inout) :: nirred_comb,iirred_comb
3144  logical,optional,intent(in) :: compute,anharmstr,spcoupling,disp
3145  integer,optional,intent(in) :: nbody
3146  logical,optional,intent(in) :: only_odd_power,only_even_power
3147 !arrays
3148  integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff_sym+nstr,ncoeff_sym+nstr)
3149  integer,intent(in) :: list_coeff(6,ncoeff_sym,nsym),list_str(nstr,nsym,2)
3150  integer,intent(in) :: index_coeff_in(power_disp-1)
3151  integer,intent(out) :: list_combination(power_disp_max,nirred_comb)
3152  character(len=5),intent(in) :: symbols(natom)
3153 !Local variables ---------------------------------------
3154 !scalar
3155  integer :: icoeff1,icoeff2,nbody_in,ii,jj,nbody_count
3156  integer :: ndisp_out,nstrain
3157  logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling
3158  logical :: need_only_odd_power,need_only_even_power,compute_sym,need_disp
3159 !arrays
3160  integer :: powers(power_disp)
3161  integer,allocatable :: index_coeff(:)
3162 ! *************************************************************************
3163 
3164  !Set the inputs
3165  need_compute = .TRUE.
3166  need_anharmstr = .TRUE.
3167  need_spcoupling = .TRUE.
3168  need_disp = .TRUE.
3169  need_only_odd_power = .FALSE.
3170  need_only_even_power = .FALSE.
3171  compute_sym = .FALSE. !Never compute the symmetric combinations here
3172  nbody_in = 0 !all kind of terms
3173  if(present(compute)) need_compute = compute
3174  if(present(nbody)) nbody_in = nbody
3175  if(present(anharmstr)) need_anharmstr = anharmstr
3176  if(present(spcoupling)) need_spcoupling = spcoupling
3177  if(present(disp)) need_disp = disp
3178  if(present(only_odd_power)) need_only_odd_power = only_odd_power
3179  if(present(only_even_power)) need_only_even_power = only_even_power
3180 
3181  if(power_disp <= power_disp_max)then
3182 
3183    !  Initialisation of variables
3184    ABI_MALLOC(index_coeff,(power_disp))
3185    index_coeff(1:power_disp-1) = index_coeff_in(:)
3186    !   write(std_out,*) "DEBUG index_coff in recursive CCL: ", index_coeff
3187    !  Loop over ncoeff+nstr
3188    do icoeff1=icoeff,ncoeff+nstr
3189 
3190      !    Reset the flag compatible and possible
3191      compatible = .TRUE.
3192      possible   = .TRUE.
3193 
3194      !    If the power_disp is one, we need to set icoeff to icoeff1
3195      if(power_disp==1) then
3196        if(icoeff1<=ncoeff .and. compatibleCoeffs(icoeff,icoeff1)==0)then
3197          ! is_displacement and compatible
3198          compatible = .FALSE.
3199        end if
3200      end if
3201      !    If the distance between the 2 coefficients is superior than the cut-off, we cycle.
3202      do icoeff2=1,power_disp-1
3203        !      write(std_out,*) "icoeff1: ", icoeff1
3204        !      write(std_out,*) "icoeff2: ", icoeff2, "index_icoeff2: ", index_coeff(icoeff2)
3205        if(icoeff1 <= ncoeff .and. index_coeff(icoeff2) <=ncoeff)then
3206          if(compatibleCoeffs(index_coeff(icoeff2),icoeff1)==0)then
3207            compatible = .FALSE.
3208          end if
3209        endif
3210      end do
3211 
3212      if (.not.compatible) cycle !The distance is not compatible
3213 
3214      !    Set the index of the new coeff in the list
3215      index_coeff(power_disp) = icoeff1
3216      !    Do some checks
3217      !    -------------
3218      !    1-Check if the coefficient is full anharmonic strain and if we need to compute it
3219      if(all(index_coeff > ncoeff))then
3220        compatible = (need_anharmstr .or. need_spcoupling)
3221        possible = need_anharmstr
3222      end if
3223      !    2-Check if the coefficient is strain-coupling and if we need to compute it
3224      if(any(index_coeff <= ncoeff) .and. any(index_coeff > ncoeff))then
3225        !       write(std_out,*) "index_coeff", index_coeff,"need_spcoupling",need_spcoupling
3226        possible   = need_spcoupling
3227        compatible = need_spcoupling
3228        if(count(index_coeff > ncoeff) > max_power_strain)then
3229          possible = .false.
3230          compatible = .false.
3231        end if
3232      end if
3233      !    3-Check if the coefficient is only disp and if we need to compute it
3234      if(all(index_coeff <= ncoeff))then
3235        !       write(std_out,*) "index_coeff", index_coeff,"need_dis",need_disp
3236        compatible = (need_disp .or. need_spcoupling)
3237        possible = need_disp
3238      end if
3239      !    4-Count number of Strain and number of displacements for compute symmetric terms
3240      nstrain = 0
3241      ndisp_out = 0
3242      do ii=1,power_disp
3243        if(index_coeff(ii) > 0 .and. index_coeff(ii) <= ncoeff)then
3244          ndisp_out = ndisp_out + 1
3245        else
3246          nstrain = nstrain +1
3247          index_coeff(ii) = index_coeff(ii) - ncoeff + ncoeff_sym
3248        end if
3249      end do
3250 
3251      if(power_disp >= power_disp_min) then
3252 
3253        !      count the number of body
3254        powers(:) = 1
3255        do ii=1,power_disp
3256          do jj=ii+1,power_disp
3257            if (powers(jj) == 0) cycle
3258            if(index_coeff(ii)==index_coeff(jj))then
3259              powers(ii) = powers(ii) + 1
3260              powers(jj) = 0
3261            end if
3262          end do
3263        end do
3264        nbody_count = count(powers /= 0)
3265        !       write(std_out,*) "powers: ", powers
3266        !       write(std_out,*) "nbody_count: ", nbody_count
3267        !      check the only_odd and only_even flags
3268        if(any(mod(powers(1:power_disp),2) /=0) .and. need_only_even_power) then
3269          possible = .false.
3270        end if
3271        if(any(mod(powers(1:power_disp),2) ==0) .and. need_only_odd_power)then
3272          possible = .false.
3273        end if
3274        !      Check the nbody flag
3275        if(nbody_in /= 0)then
3276          if(power_disp-count(powers==0) > nbody_in) then
3277            possible = .false.
3278            compatible = .false.
3279          end if
3280        end if
3281 
3282        if(possible) then
3283          !        increase coefficients and set it
3284          !         nmodel_tot = nmodel_tot + 1
3285          !         if(need_compute)then
3286          !           list_combination(1:power_disp,nmodel_tot) = index_coeff
3287          !         end if
3288          !nmodel_tot_test = 0
3289          !         !Start from second symmetry in Symmetric Combinations
3290          !         isym_in_test = 2
3291          !         idisp_in_test = power_disp
3292          !         ndisp_test = power_disp
3293          !         index_coeff_tmp = index_coeff
3294          !Count anharmonic strain terms
3295          if(ndisp_out == 0 .and. nstrain > 0)then
3296            nirred_comb = nirred_comb +1
3297            iirred_comb = iirred_comb +1
3298            if(need_compute)then
3299              list_combination(1:power_disp,iirred_comb) = index_coeff
3300            endif
3301          else !Else counst symmetric terms of atomic displacement (pure disp or disp/strain)
3302            !Store index for each combination of irreducible terms to later parallely compute symmetric combinations
3303            nirred_comb = nirred_comb +1
3304            iirred_comb = iirred_comb +1
3305            !write(std_out,*) "DEBUG index_coeff: ", index_coeff
3306            !write(std_out,*) "DEBUG ndisp_out: ", ndisp_out
3307            !write(std_out,*) "DEBUG nstrain: ", nstrain
3308            !write(std_out,*) "DEBUG nmodel_start: ", nmodel_start
3309            if(need_compute)then
3310              list_combination(:,iirred_comb) = 0
3311              list_combination(:ndisp_out+nstrain,iirred_comb) = index_coeff
3312            endif
3313          end if !ndisp_out == 0 .and.n nstrain >0
3314        end if!possible
3315      end if!end if power_disp < power_disp_min
3316 
3317      !Change back to irreducible terms ncoeff_limit
3318      do ii=1,power_disp
3319        if(index_coeff(ii) > ncoeff_sym)then
3320          index_coeff(ii) = index_coeff(ii) + ncoeff - ncoeff_sym
3321        end if
3322      end do
3323 
3324      !    If the model is still compatbile with the input flags, we continue.
3325      if(compatible)then
3326        call computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
3327          &                                     index_coeff,list_combination,icoeff1,max_power_strain,&
3328          &                                     natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp+1,&
3329          &                                     power_disp_min,power_disp_max,symbols,comm,nbody=nbody_in,&
3330          &                                     compute=need_compute,anharmstr=need_anharmstr,&
3331          &                                     spcoupling=need_spcoupling,only_odd_power=need_only_odd_power,&
3332          &                                     only_even_power=need_only_even_power,disp=need_disp)
3333      end if
3334    end do
3335    ABI_FREE(index_coeff)
3336  end if
3337 
3338 end subroutine computeCombinationFromList

m_polynomial_coeff/computeNorder [ Functions ]

[ Top ] [ Functions ]

NAME

 computeNorder

FUNCTION

 Recursive routine to compute the order N of a all the possible coefficient
 from the list list_symcoeff and list_symstr.

INPUTS

 cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
 compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1
 list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 list_symstr(nstr_sym,nsym) = array with the list of the strain  and the symmetrics
 index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0)
 icoeff = current indexes of the cofficients (start we 1)
 icoeff_tot = current number of coefficients already computed (start we 0)
 natom = number of atoms in the unit cell
 nstr = number of coefficient for related to the atomic displacment into list_symcoeff
 nstr = number of coefficient for related to the strain into list_symstr
 ncoeff_out = number of maximum coefficients
 nrpt = number of cell
 nsym = number of symmetries in the system
 power_disp = initial power_disp to be computed (can be < power_disp_min,
              this routine will skip the firts power_disp)
 power_disp_min = minimal power_disp to be computed
 power_disp_max = maximum power_disp to be computed
 symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
 nbody = optional, number of body for the coefficients, for example:
                   0 => all the terms
                   1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp  ...
 compute = logical, optional: TRUE if we store the coefficients
                              FALSE just to count the number of coefficient
 anharmstr = logical, optional : TRUE, the anharmonic strain are computed
                                   FALSE, (default) the anharmonic strain are not computed
 distributed = logical, optional : True, the coefficients will be distributed on the CPU

OUTPUT

 icoeff = current indexes of the cofficients (start we 1)
 icoeff_tot = current number of coefficients already computed (start we 0)
 polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
                                                              the polynomial_coeff

SOURCE

2883 recursive subroutine computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,&
2884 &                                  index_coeff_in,icoeff,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,&
2885 &                                  nrpt,nsym,power_disp,power_disp_min,power_disp_max,symbols,nbody,&
2886 &                                  compute,anharmstr,spcoupling,distributed)
2887 
2888  implicit none
2889 
2890 !Arguments ---------------------------------------------
2891 !scalar
2892  integer,intent(in) :: natom,ncoeff,power_disp,power_disp_min,power_disp_max,ncoeff_out,nsym,nrpt,nstr,icoeff
2893  integer,intent(inout) :: icoeff_tot
2894  logical,optional,intent(in) :: compute,anharmstr,spcoupling,distributed
2895  integer,optional,intent(in) :: nbody
2896 !arrays
2897  integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff+nstr,ncoeff+nstr)
2898  integer,intent(in) :: list_coeff(6,ncoeff,nsym),list_str(nstr,nsym,2)
2899  integer,intent(in) :: index_coeff_in(power_disp-1)
2900  type(polynomial_coeff_type),intent(inout) :: coeffs_out(ncoeff_out)
2901  character(len=5),intent(in) :: symbols(natom)
2902 !Local variables ---------------------------------------
2903 !scalar
2904  integer :: ia,ib,ii,icoeff1,icoeff_tmp
2905  integer :: iterm,nbody_in,ncoeff_max,pa,pb
2906  integer :: ndisp_max,nterm_max
2907  real(dp):: coefficient
2908  logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling,need_distributed
2909 !arrays
2910  integer,allocatable :: index_coeff(:)
2911  character(len=200):: name
2912  type(polynomial_term_type),dimension(:),allocatable :: terms
2913  type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
2914 ! *************************************************************************
2915 
2916 !Set the inputs
2917  need_compute = .TRUE.
2918  need_anharmstr = .TRUE.
2919  need_spcoupling = .TRUE.
2920  need_distributed = .FALSE.
2921  nbody_in = 0 !all kind of terms
2922  if(present(compute)) need_compute = compute
2923  if(present(nbody)) nbody_in = nbody
2924  if(present(anharmstr)) need_anharmstr = anharmstr
2925  if(present(spcoupling)) need_spcoupling = spcoupling
2926  if(present(distributed)) need_distributed  = distributed
2927  if(power_disp <= power_disp_max)then
2928 
2929 !  Initialisation of variables
2930    nterm_max  = nsym
2931    ncoeff_max = (ncoeff+nstr)
2932    ndisp_max = power_disp
2933    icoeff_tmp = 0
2934    ABI_MALLOC(coeffs_tmp,(ncoeff_max))
2935    ABI_MALLOC(terms,(nterm_max))
2936    ABI_MALLOC(index_coeff,(power_disp))
2937 
2938    index_coeff(1:power_disp-1) = index_coeff_in(:)
2939 
2940    do icoeff1=icoeff,ncoeff+nstr
2941 !    If the distance between the 2 coefficients is superior than the cut-off,
2942 !    we cycle
2943 !    If the power_disp is one check if icoeff1 is compatible with itself
2944      if(power_disp==1) then
2945        if(icoeff1 <= ncoeff .and. compatibleCoeffs(icoeff1,icoeff1)==0)then
2946          cycle
2947        end if
2948      end if
2949      if(compatibleCoeffs(icoeff,icoeff1)==0) cycle
2950 
2951 !    Reset the flag compatible and possible
2952      compatible = .TRUE.
2953      possible   = .TRUE.
2954 
2955      index_coeff(power_disp) = icoeff1
2956      iterm = 0
2957      coefficient = one
2958 
2959      if(power_disp >= power_disp_min) then
2960        block
2961          logical :: reverse(ndisp_max)
2962          reverse(:) = .False.
2963          call generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,&
2964            &                                 ndisp_max,nrpt,nstr,nsym,iterm,terms, reverse=reverse)
2965        end block
2966 
2967        if(iterm > 0)then
2968 !        Do some checks
2969 !        -------------
2970 !        1-Check if the coefficient is full anharmonic strain and if we need to compute it
2971          if(terms(1)%ndisp == 0)then
2972            compatible = (need_anharmstr .or. need_spcoupling)
2973            possible = need_anharmstr
2974          end if
2975 !        1-Check if the coefficient is strain-coupling and if we need to compute it
2976          if(terms(1)%nstrain > 0.and.terms(1)%ndisp > 0)then
2977            possible   = need_spcoupling
2978            compatible = need_spcoupling
2979          end if
2980 
2981 !        ------------
2982 !        2-Check if this terms is compatible with nbody
2983          if(nbody_in > 0)then
2984            pa = 1 ; pb = 1
2985            ia = 0 ; ib = 0
2986 !          Count the number of terms and the power_disp
2987            do ii=1,terms(1)%ndisp
2988              if(terms(1)%nstrain > 0) then
2989                pb = pb*terms(1)%power_disp(ii)
2990                ib = ib + 1
2991              else
2992                pa = pa*terms(1)%power_disp(ii)
2993                ia = ia + 1
2994              end if
2995            end do
2996            if(ia <= nbody_in)then
2997              if(ia==nbody_in.and.abs(mod(pa,2)) < tol16)then
2998                if(ib==0)then
2999                  compatible = .FALSE.
3000                  possible   = .TRUE.
3001                else if (ib==nbody_in.and.abs(mod(pb,2)) < tol16) then
3002                  compatible = .FALSE.
3003                  possible   = .TRUE.
3004                else
3005                 possible = .FALSE.
3006                 compatible = .FALSE.
3007                end if
3008              else
3009                 possible = .FALSE.
3010                 compatible = .FALSE.
3011              end if
3012            else
3013              compatible = .FALSE.
3014              possible = .FALSE.
3015            end if
3016          end if
3017 
3018          if(possible)then
3019 !          increase coefficients and set it
3020            icoeff_tmp = icoeff_tmp + 1
3021            icoeff_tot = icoeff_tot + 1
3022            call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),&
3023 &                                     terms(1:iterm),check=.true.)
3024 
3025          end if
3026        end if
3027 
3028 !      Deallocate the terms
3029        do iterm=1,nterm_max
3030          call polynomial_term_free(terms(iterm))
3031        end do
3032 
3033      end if!end if power_disp < power_disp_min
3034 
3035      if(compatible)then
3036        call computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,index_coeff,&
3037 &                         icoeff1,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,nrpt,nsym,power_disp+1,&
3038 &                         power_disp_min,power_disp_max,symbols,nbody=nbody_in,compute=need_compute,&
3039 &                         anharmstr=need_anharmstr,spcoupling=need_spcoupling)
3040      end if
3041    end do
3042 
3043    ABI_FREE(terms)
3044    ABI_FREE(index_coeff)
3045 
3046 !  Transfer in the final array
3047    icoeff1 = 0
3048    do icoeff_tmp=1,ncoeff_max
3049      if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3050 !      Increase icoeff and fill the coeffs_out array
3051        icoeff_tot = icoeff_tot + 1
3052        if(need_compute)then
3053          name = ''
3054 !        Get the name of this coefficient
3055          call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
3056          call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
3057 &                                   coeffs_out(icoeff_tot),coeffs_tmp(icoeff_tmp)%terms,&
3058 &                                   name=name)
3059        end if
3060      end if
3061    end do
3062 !  Deallocation
3063    do icoeff1=1,ncoeff_max
3064      call polynomial_coeff_free(coeffs_tmp(icoeff1))
3065    end do
3066    ABI_FREE(coeffs_tmp)
3067  end if
3068 
3069 end subroutine computeNorder

m_polynomial_coeff/generateTermsFromList [ Functions ]

[ Top ] [ Functions ]

NAME

 generateTermsFromList

FUNCTION

 Compute for a given list of index the correspondig set of terms

INPUTS

 cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
 index_coeff_in(ndisp) = list of coefficients to be computed
 list_symcoeff(6,ncoeff,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 list_symstr(nstr,nsym) = array with the list of the strain  and the symmetrics
 ncoeff = number of maximum coefficients in the list_symcoeff
 ndisp = number of maximum diplacement (phonon + strain)
 nrpt = number of cell
 nsym = number of symmetries in the system

OUTPUT

 terms<(type(polynomial_term_type)>(nterm)  = list of terms
 nterm = number of ouput terms

SOURCE

3645 subroutine generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,ndisp_max,&
3646 &                                nrpt,nstr,nsym,nterm,terms, reverse)
3647 
3648  implicit none
3649 
3650 !Arguments ------------------------------------
3651 !scalar
3652  integer,intent(in) :: ndisp_max,ncoeff,nrpt,nstr,nsym
3653  integer,intent(out):: nterm
3654 !arrays
3655  integer,intent(in) :: index_coeff(ndisp_max)
3656  integer,intent(in) :: cell(3,nrpt),list_coeff(6,ncoeff,nsym)
3657  integer,intent(in) :: list_str(nstr,nsym,2)
3658  logical, intent(in) :: reverse(ndisp_max)
3659  type(polynomial_term_type),intent(out) :: terms(nsym)
3660 !Local variables-------------------------------
3661 !scalar
3662  integer :: ia,ib,icoeff_str,idisp,irpt
3663  integer :: isym,ndisp,nstrain,mu
3664  real(dp):: weight
3665 !arrays
3666  integer :: atindx(2,ndisp_max),cells(3,2,ndisp_max),dir_int(ndisp_max),strain(ndisp_max)
3667  integer :: power_disps(ndisp_max),power_strain(ndisp_max)
3668 
3669 ! *************************************************************************
3670  nterm = 0
3671 !Loop over symetries
3672  do isym=1,nsym
3673 !Treat this coeff
3674    weight = 1
3675    ndisp = 0
3676    nstrain = 0
3677    do idisp=1,ndisp_max
3678 !    Get index of this displacement term
3679 !    Check if the index is not zero
3680      if(index_coeff(idisp)==0)cycle
3681      if(index_coeff(idisp)<=ncoeff)then
3682        ndisp = ndisp + 1
3683        mu   = list_coeff(1,index_coeff(idisp),isym)
3684        if( reverse(idisp)) then
3685          ia   = list_coeff(3,index_coeff(idisp),isym)
3686          ib   = list_coeff(2,index_coeff(idisp),isym)
3687          irpt = list_coeff(4,index_coeff(idisp),isym)
3688          irpt = find_opposite_irpt(cell, irpt)
3689          weight = -weight*list_coeff(5,index_coeff(idisp),isym)
3690        else
3691          ia   = list_coeff(2,index_coeff(idisp),isym)
3692          ib   = list_coeff(3,index_coeff(idisp),isym)
3693          irpt = list_coeff(4,index_coeff(idisp),isym)
3694          weight = weight*list_coeff(5,index_coeff(idisp),isym)
3695        end if
3696 !      Fill First term arrays
3697        atindx(1,idisp) = ia; atindx(2,idisp) = ib;
3698        dir_int(idisp) = mu
3699        power_disps(idisp)   = 1
3700        cells(:,1,idisp) = (/0,0,0/)
3701        cells(:,2,idisp) = cell(:,irpt)
3702      else
3703        nstrain = nstrain + 1
3704        icoeff_str = index_coeff(idisp)-ncoeff
3705        strain(nstrain) = list_str(icoeff_str,isym,1)
3706        power_strain(nstrain)  = 1
3707        weight = weight*list_str(icoeff_str,isym,2)
3708      end if
3709    end do
3710    nterm = nterm + 1
3711    call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(nterm),power_disps,&
3712 &                            power_strain,strain,weight,check=.true.)
3713  end do!end do sym
3714 
3715 
3716 end subroutine generateTermsFromList

m_polynomial_coeff/getCoeffFromList [ Functions ]

[ Top ] [ Functions ]

NAME

 getCoeffFromList

FUNCTION

 get the index of a coefficient into the list_coeff

INPUTS

 list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 ia = index of the atom 1
 ib = index of the atom 1
 irpt = indexes of the cell of the second atom
 mu = direction of the IFC
 ncoeff = number of total coefficients in the list

OUTPUT

 coeff = index of the coefficient

SOURCE

3579 function getCoeffFromList(list_coeff,ia,ib,irpt,mu,ncoeff) result(coeff)
3580 
3581  implicit none
3582 
3583 !Arguments ------------------------------------
3584 !scalar
3585  integer,intent(in) :: ia,ib,irpt,mu,ncoeff
3586  integer :: coeff
3587 !arrays
3588  integer,intent(in) :: list_coeff(6,ncoeff)
3589 !Local variables-------------------------------
3590 !scalar
3591  integer :: icoeff
3592 !arrays
3593 
3594 ! *************************************************************************
3595  coeff = 0
3596  !write(std_out,*) "DEBUG shape inside FromList: ", shape(list_coeff)
3597  do icoeff = 1,ncoeff
3598    if(mu==list_coeff(1,icoeff).and.&
3599 &     ia==list_coeff(2,icoeff).and.&
3600 &     ib==list_coeff(3,icoeff).and.&
3601 &     irpt==list_coeff(4,icoeff))then!.and.&
3602 !&     abs(weight-list_coeff(5,icoeff)) < tol16) then
3603      coeff = icoeff
3604      exit
3605    end if
3606  end do
3607 
3608 end function getCoeffFromList

m_polynomial_coeff/polynomial_coeff_broadcast [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_broadcast

FUNCTION

  MPI broadcast  polynomial_coefficent datatype

INPUTS

  source = rank of source
  comm = MPI communicator

SIDE EFFECTS

  coefficients<type(polynomial_coefficent_type)>= Input if node is source,
                              other nodes returns with a completely initialized instance.

SOURCE

606 subroutine polynomial_coeff_broadcast(coefficients, source, comm)
607 
608  implicit none
609 
610 !Arguments ------------------------------------
611 !array
612  type(polynomial_coeff_type),intent(inout) :: coefficients
613  integer, intent(in) :: source,comm
614 
615 !Local variables-------------------------------
616 !scalars
617  integer :: ierr,ii
618 !arrays
619 
620 ! *************************************************************************
621 
622 
623  if (xmpi_comm_size(comm) == 1) return
624 
625 ! Free the output
626  if (xmpi_comm_rank(comm) /= source) then
627    call polynomial_coeff_free(coefficients)
628  end if
629 
630  ! Transmit variables
631   call xmpi_bcast(coefficients%name, source, comm, ierr)
632   call xmpi_bcast(coefficients%nterm, source, comm, ierr)
633   call xmpi_bcast(coefficients%coefficient, source, comm, ierr)
634 
635  !Allocate arrays on the other nodes.
636   if (xmpi_comm_rank(comm) /= source) then
637     ABI_MALLOC(coefficients%terms,(coefficients%nterm))
638     do ii=1,coefficients%nterm
639       call polynomial_term_free(coefficients%terms(ii))
640     end do
641   end if
642 ! Set the number of term on each node (needed for allocations of array)
643   do ii = 1,coefficients%nterm
644     call xmpi_bcast(coefficients%terms(ii)%ndisp, source, comm, ierr)
645     call xmpi_bcast(coefficients%terms(ii)%nstrain, source, comm, ierr)
646   end do
647 
648 ! Allocate arrays on the other nodes
649   if (xmpi_comm_rank(comm) /= source) then
650     do ii = 1,coefficients%nterm
651       ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
652       coefficients%terms(ii)%atindx = 0
653       ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
654       ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
655       ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
656       ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
657       ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
658     end do
659   end if
660 
661 ! Transfert value
662   do ii = 1,coefficients%nterm
663       call xmpi_bcast(coefficients%terms(ii)%weight, source, comm, ierr)
664       call xmpi_bcast(coefficients%terms(ii)%atindx, source, comm, ierr)
665       call xmpi_bcast(coefficients%terms(ii)%direction, source, comm, ierr)
666       call xmpi_bcast(coefficients%terms(ii)%cell, source, comm, ierr)
667       call xmpi_bcast(coefficients%terms(ii)%power_disp, source, comm, ierr)
668       call xmpi_bcast(coefficients%terms(ii)%power_strain, source, comm, ierr)
669       call xmpi_bcast(coefficients%terms(ii)%strain, source, comm, ierr)
670   end do
671 
672 
673 end subroutine polynomial_coeff_broadcast

m_polynomial_coeff/polynomial_coeff_evaluate [ Functions ]

[ Top ] [ Functions ]

NAME

  polynomial_coeff_evaluate

FUNCTION

  Compute the energy related to the coefficients from
  fitted polynome

INPUTS

  coefficients(ncoeff)<type(polynomial_coeff_type)> = list of coefficients
  disp(3,natom_sc) = atomics displacement between configuration and the reference
  natom_sc = number of atoms in the supercell
  natom_uc = number of atoms in the unit cell
  ncoeff   = number of coefficients
  sc_size(3) = size of the supercell (2 2 2 for example)
  strain(6) = strain between configuration and the reference
  cells(ncell) = number of the cells into the supercell (1,2,3,4,5)
  ncell   = total number of cell to treat by this cpu
  index_cells(3,ncell) = indexes of the cells into  supercell (-1 -1 -1 ,...,1 1 1)
  comm=MPI communicator

OUTPUT

  energy = contribution to the energy
  energy_coeff(ncoeff) = energy contribution of each anharmonic term
  fcart(3,natom) = contribution  to the forces
  strten(6) = contribution to the stress tensor

SOURCE

1018 subroutine polynomial_coeff_evaluate(coefficients,disp,energy,energy_coeff,fcart,natom_sc,natom_uc,ncoeff,sc_size,&
1019 &                                    strain,strten,ncell,index_cells,comm,filename)
1020 
1021 !Arguments ------------------------------------
1022 ! scalar
1023   real(dp),intent(out):: energy
1024   integer, intent(in) :: ncell,ncoeff,natom_sc,natom_uc
1025   integer, intent(in) :: comm
1026   character(len=fnlen),optional,intent(in) :: filename
1027 ! array
1028   real(dp),intent(out):: strten(6)
1029   real(dp),intent(in) :: strain(6)
1030   real(dp),intent(out):: fcart(3,natom_sc)
1031   real(dp),intent(in) :: disp(3,natom_sc)
1032   real(dp),optional,intent(out):: energy_coeff(ncoeff)
1033   integer,intent(in) :: index_cells(4,ncell)
1034   integer,intent(in) :: sc_size(3)
1035   type(polynomial_coeff_type),intent(in) :: coefficients(ncoeff)
1036  !Local variables-------------------------------
1037 ! scalar
1038   integer :: i1,i2,i3,ia1,ib1,ia2,ib2,idir1,idir2,ierr,ii
1039   integer :: icoeff,iterm,idisp1,idisp2,idisp1_strain,idisp2_strain,icell,ndisp
1040   integer :: nstrain,ndisp_tot,power_disp,power_strain,unit_out
1041   real(dp):: coeff,disp1,disp2,tmp1,tmp2,tmp3,weight
1042   logical :: file_opened
1043 ! array
1044   integer :: cell_atoma1(3),cell_atoma2(3)
1045   integer :: cell_atomb1(3),cell_atomb2(3)
1046   character(len=500) :: msg
1047   character(len=fnlen) :: name_file
1048 ! *************************************************************************
1049 
1050 ! Check
1051   if (any(sc_size <= 0)) then
1052     write(msg,'(a,a)')' No supercell found for getEnergy'
1053     ABI_ERROR(msg)
1054   end if
1055 
1056   if(present(filename)) name_file = filename
1057 
1058 ! Initialisation of variables
1059   energy     = zero
1060   fcart(:,:) = zero
1061   strten(:)  = zero
1062   energy_coeff(:) = zero
1063   do icell = 1,ncell
1064     ii = index_cells(4,icell);
1065     i1=index_cells(1,icell); i2=index_cells(2,icell); i3=index_cells(3,icell)
1066     ia1 = 0 ; ib1 = 0
1067 !   Loop over coefficients
1068     do icoeff=1,ncoeff
1069 !     Set the value of the coefficient
1070       coeff = coefficients(icoeff)%coefficient
1071 !     Loop over terms of this coefficient
1072       do iterm=1,coefficients(icoeff)%nterm
1073 !       Set the weight of this term
1074         weight =coefficients(icoeff)%terms(iterm)%weight
1075         tmp1 = one
1076         ndisp = coefficients(icoeff)%terms(iterm)%ndisp
1077         nstrain = coefficients(icoeff)%terms(iterm)%nstrain
1078         ndisp_tot = ndisp + nstrain
1079 !       Loop over displacement and strain
1080         do idisp1=1,ndisp_tot
1081 !         Set to one the acculation of forces and strain
1082           tmp2 = one
1083           tmp3 = zero
1084 
1085 !         Strain case idisp > ndisp
1086           if (idisp1 > ndisp)then
1087             tmp3 = one
1088 !           Set the power_strain of the strain:
1089             idisp1_strain = idisp1 - ndisp
1090             power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain)
1091 !           Get the direction of the displacement or strain
1092             idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain)
1093             if(abs(strain(idir1)) > tol10)then
1094 !             Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
1095               tmp1 = tmp1 * (strain(idir1))**power_strain
1096               if(power_strain > 1) then
1097 !               Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
1098                 tmp3 = tmp3 *  power_strain*(strain(idir1))**(power_strain-1)
1099               end if
1100             else
1101               tmp1 = zero
1102               if(power_strain > 1) then
1103                 tmp3 = zero
1104               end if
1105             end if
1106           else
1107 !           Set the power_disp of the displacement:
1108             power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1)
1109 !           Get the direction of the displacement or strain
1110             idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1)
1111 !           Displacement case idir = 1, 2  or 3
1112 !           indexes of the cell of the atom a
1113             cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1)
1114             if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then
1115 !             if the cell is not 0 0 0 we apply PBC:
1116               cell_atoma1(1) =  i1 + cell_atoma1(1)
1117               cell_atoma1(2) =  i2 + cell_atoma1(2)
1118               cell_atoma1(3) =  i3 + cell_atoma1(3)
1119               call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3))
1120 !             index of the first atom (position in the supercell if the cell is not 0 0 0)
1121               ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1122 &                   (cell_atoma1(2)-1)*sc_size(3)*natom_uc+&
1123 &                   (cell_atoma1(3)-1)*natom_uc+&
1124 &                   coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
1125             else
1126 !             index of the first atom (position in the supercell if the cell is 0 0 0)
1127               ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
1128             end if
1129 
1130 !           indexes of the cell of the atom b  (with PBC) same as ia1
1131             cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1)
1132             if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then
1133               cell_atomb1(1) =  i1 + cell_atomb1(1)
1134               cell_atomb1(2) =  i2 + cell_atomb1(2)
1135               cell_atomb1(3) =  i3 + cell_atomb1(3)
1136               call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3))
1137 
1138 !             index of the second atom in the (position in the supercell  if the cell is not 0 0 0)
1139               ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1140 &                   (cell_atomb1(2)-1)*sc_size(3)*natom_uc+&
1141 &                   (cell_atomb1(3)-1)*natom_uc+&
1142 &                   coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
1143             else
1144 !             index of the first atom (position in the supercell if the cell is 0 0 0)
1145               ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
1146             end if
1147 
1148 !           Get the displacement for the both atoms
1149             disp1 = disp(idir1,ia1)
1150             disp2 = disp(idir1,ib1)
1151 
1152             if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then
1153 !           Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
1154               tmp1 = tmp1 * (disp1-disp2)**power_disp
1155               if(power_disp > 1) then
1156 !               Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...))
1157                 tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1)
1158               end if
1159             else
1160               tmp1 = zero
1161               if(power_disp > 1) then
1162                 tmp2 = zero
1163               end if
1164             end if
1165           end if
1166 
1167           do idisp2=1,ndisp_tot
1168 
1169             if(idisp2 /= idisp1) then
1170               if (idisp2 > ndisp)then
1171                 idisp2_strain = idisp2 - ndisp
1172                 idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain)
1173 !               Strain case
1174 !               Set the power_strain of the strain:
1175                 power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain)
1176 !               Accumulate energy forces
1177                 tmp2 = tmp2 * (strain(idir2))**power_strain
1178 !               Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
1179                 tmp3 = tmp3 * (strain(idir2))**power_strain
1180               else
1181                 idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2)
1182                 cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2)
1183                 if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then
1184                   cell_atoma2(1) =  i1 + cell_atoma2(1)
1185                   cell_atoma2(2) =  i2 + cell_atoma2(2)
1186                   cell_atoma2(3) =  i3 + cell_atoma2(3)
1187                   call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3))
1188 !                 index of the first atom (position in the supercell and direction)
1189 !                 if the cell of the atom a is not 0 0 0 (may happen)
1190                   ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1191 &                       (cell_atoma2(2)-1)*sc_size(3)*natom_uc+&
1192 &                       (cell_atoma2(3)-1)*natom_uc+&
1193 &                       coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
1194                 else
1195 !                 index of the first atom (position in the supercell and direction)
1196                   ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
1197                 end if
1198 
1199                 cell_atomb2= coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2)
1200 
1201                 if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then
1202 !                 indexes of the cell2 (with PBC)
1203                   cell_atomb2(1) =  i1 + cell_atomb2(1)
1204                   cell_atomb2(2) =  i2 + cell_atomb2(2)
1205                   cell_atomb2(3) =  i3 + cell_atomb2(3)
1206                   call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3))
1207 
1208 !                 index of the second atom in the (position in the supercell)
1209                   ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
1210 &                       (cell_atomb2(2)-1)*sc_size(3)*natom_uc+&
1211 &                       (cell_atomb2(3)-1)*natom_uc+&
1212 &                       coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
1213                 else
1214                   ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
1215                 end if
1216 
1217                 disp1 = disp(idir2,ia2)
1218                 disp2 = disp(idir2,ib2)
1219 
1220 !               Set the power_disp of the displacement:
1221                 power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2)
1222                 tmp2 = tmp2 * (disp1-disp2)**power_disp
1223                 tmp3 = tmp3 * (disp1-disp2)**power_disp
1224 
1225               end if
1226             end if
1227           end do
1228 
1229           if(idisp1 > ndisp)then
1230 !           Accumule stress tensor
1231             strten(idir1) = strten(idir1) + coeff * weight * tmp3
1232           else
1233 !           Accumule  forces
1234             fcart(idir1,ia1) =  fcart(idir1,ia1)  + coeff * weight * tmp2
1235             fcart(idir1,ib1) =  fcart(idir1,ib1)  - coeff * weight * tmp2
1236           end if
1237         end do
1238 
1239         energy_coeff(icoeff) = energy_coeff(icoeff) + coeff * weight * tmp1
1240 !       accumule energy
1241         energy = energy +  coeff * weight * tmp1
1242 
1243       end do
1244     end do
1245   end do
1246 
1247 
1248 ! MPI_SUM
1249   call xmpi_sum(energy, comm, ierr)
1250   call xmpi_sum(fcart , comm, ierr)
1251   call xmpi_sum(strten , comm, ierr)
1252 
1253 !Write to anharmonic_energy_terms.out ORIGINAL
1254   INQUIRE(FILE=name_file,OPENED=file_opened,number=unit_out)
1255   if(file_opened .eqv. .TRUE.)then
1256     do icoeff=1,ncoeff
1257       call xmpi_sum(energy_coeff(icoeff), comm, ierr)
1258      ! Marcus write energy contributions of anharmonic terms to file
1259       if(icoeff <ncoeff)then
1260         write(unit_out,'(A,1ES24.16)',advance='no')  '    ',energy_coeff(icoeff)
1261       else if(icoeff==ncoeff)then
1262         write(unit_out,'(A,1ES24.16)',advance='yes') '    ',energy_coeff(icoeff)
1263       end if
1264     enddo
1265   end if
1266 
1267 
1268 end subroutine polynomial_coeff_evaluate

m_polynomial_coeff/polynomial_coeff_free [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_free

FUNCTION

 Free polynomial_coeff datatype

INPUTS

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

OUTPUT

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

SOURCE

308 subroutine polynomial_coeff_free(polynomial_coeff)
309 
310  implicit none
311 
312 !Arguments ------------------------------------
313 !scalars
314 !arrays
315  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
316 !Local variables-------------------------------
317 !scalar
318  integer :: ii
319 !arrays
320 
321 ! *************************************************************************
322 
323  if(allocated(polynomial_coeff%terms))then
324    do ii = 1,polynomial_coeff%nterm
325      call polynomial_term_free(polynomial_coeff%terms(ii))
326    end do
327    ABI_FREE(polynomial_coeff%terms)
328  end if
329  polynomial_coeff%name = ""
330  polynomial_coeff%nterm = 0
331  polynomial_coeff%coefficient = zero
332 
333 
334 end subroutine polynomial_coeff_free

m_polynomial_coeff/polynomial_coeff_getEvenAnhaStrain [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_getEvenAnhaStrain

FUNCTION

 Get even anharmonic strain terms in defined range of order

INPUTS

OUTPUT

 polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
                                                              the polynomial_coeff
 ncoeff_out = number of coefficients

SOURCE

3934 subroutine polynomial_coeff_getEvenAnhaStrain(strain_terms,crystal,irred_ncoeff,power_strain,comm)
3935 
3936  implicit none
3937 
3938 !Arguments ------------------------------------
3939 type(polynomial_coeff_type),allocatable,intent(inout) :: strain_terms(:)
3940 type(crystal_t), intent(inout) :: crystal
3941 integer,intent(out) :: irred_ncoeff
3942 integer,intent(in) :: power_strain(2)
3943 integer,intent(in) :: comm
3944 !scalars
3945 !arrays
3946 !Local variables-------------------------------
3947 real(dp) :: cutoff,coeff_ini
3948 integer :: ncoeff
3949 integer :: power_strph
3950 integer :: option
3951 integer :: icoeff1,icoeff2,start
3952 integer :: ncoeff_out
3953 logical,allocatable :: same(:)
3954 type(polynomial_coeff_type),allocatable :: strain_terms_tmp(:)
3955 !scalar
3956 !arrays
3957 integer :: sc_size(3)
3958 ! *************************************************************************
3959 
3960 !Initialize Variables for call to polynomial_coeff_getNorder
3961 cutoff = zero
3962 power_strph = zero
3963 option = 0
3964 sc_size = (/1,1,1/)
3965 coeff_ini = 1000000
3966 
3967 
3968 call polynomial_coeff_getNorder(strain_terms_tmp,crystal,cutoff,ncoeff,ncoeff_out,power_strain,&
3969 &                               power_strph,option,sc_size,comm,anharmstr=.true.,spcoupling=.false.,&
3970 &                               only_odd_power=.false.,only_even_power=.true.,compute_symmetric=.false.,&
3971                                 verbose=.false.)
3972 
3973 
3974 !TODO Probably put in one routine
3975 !Get irreducible strain
3976 ABI_MALLOC(same,(ncoeff_out))
3977 same = .false.
3978 do icoeff1 =1,ncoeff_out
3979         start = icoeff1 + 1
3980         !write(*,*) "name coeff_",icoeff1,": ", strain_terms_tmp(icoeff1)%name
3981         do icoeff2=start,ncoeff_out
3982                 if(.not.same(icoeff2)) same(icoeff2) = coeffs_compare(strain_terms_tmp(icoeff1),strain_terms_tmp(icoeff2))
3983         enddo
3984         !write(*,*) "same(",icoeff1,"): ", same(icoeff1)
3985 enddo
3986 
3987 irred_ncoeff = ncoeff_out - count(same)
3988 ABI_MALLOC(strain_terms,(irred_ncoeff))
3989 
3990 !Transfer irreducible strain to output array
3991 icoeff2=0
3992 icoeff1=0
3993 do icoeff1=1,ncoeff_out
3994         if(.not.same(icoeff1))then
3995                 icoeff2=icoeff2 + 1
3996                 call polynomial_coeff_init(coeff_ini,strain_terms_tmp(icoeff1)%nterm,strain_terms(icoeff2),&
3997 &               strain_terms_tmp(icoeff1)%terms,strain_terms_tmp(icoeff1)%name,check=.TRUE.)
3998         endif
3999 enddo
4000 
4001 !TEST MS write coefficients to xml to check
4002 !fname="EvenAnhaStrain.xml"
4003 !call polynomial_coeff_writeXML(strain_terms,irred_ncoeff,filename=fname,newfile=.true.)
4004 !write(*,*) "Strain terms to xml done"
4005 
4006 
4007 !Deallocateion
4008 call polynomial_coeff_list_free(strain_terms_tmp)
4009 ABI_FREE(same)
4010 
4011 end subroutine polynomial_coeff_getEvenAnhaStrain

m_polynomial_coeff/polynomial_coeff_getList [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_getList

FUNCTION

 Get the list of all  the possible coefficients for the polynome

INPUTS

 cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
 dist(3,natom,natom,nrpt) = distance between atoms atm1 is in the cell 0 0 0
                                                   atm2 is in the nrpt cell (see cell(3,nrpt))
                            for each component x,y and z
 crystal<type(crystal_t)> = datatype with all the information for the crystal
 natom = number of atoms in the unit cell
 nrpt  = number of cell in the supercell

OUTPUT

 list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 list_symstr(nstr_sym,nsym) = array with the list of the strain  and the symmetrics
 nstr_sym = number of coefficient for the strain
 ncoeff_sym = number of coefficient for the IFC
 range_ifc(3) = maximum cut-off for the inter atomic forces constants in each direction
 sc_size(3) = optional,size of the supercell used for the fit.
               For example if you want to fit 2x2x2 cell the interation
               Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process
               If check_pbc is true we remove these kind of terms

SOURCE

1423 subroutine polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1424 &                                   natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size,&
1425 &                                   fit_iatom)
1426 
1427  implicit none
1428 
1429 !Arguments ------------------------------------
1430 !scalars
1431  integer,intent(in) :: natom,nrpt
1432  real(dp), intent(in)   :: cutoff
1433  integer,intent(out) :: ncoeff_sym,nstr_sym
1434  integer,optional,intent(in):: fit_iatom
1435 !arrays
1436  integer,intent(in) :: cell(3,nrpt)
1437  real(dp),intent(in):: dist(3,natom,natom,nrpt)
1438  type(crystal_t), intent(in) :: crystal
1439  integer,allocatable,intent(out) :: list_symcoeff(:,:,:),list_symstr(:,:,:)
1440  integer,optional,intent(in) :: sc_size(3)
1441  real(dp),intent(in):: range_ifc(3)
1442 !Local variables-------------------------------
1443 !scalar
1444  integer :: ia,ib,icoeff,icoeff2,icoeff_tot,icoeff_tmp,idisy1,idisy2,ii
1445  integer :: ipesy1,ipesy2,isym,irpt,irpt3,irpt_ref,irpt_sym
1446  integer :: jj,jsym,mu,fit_iatom_in
1447  integer :: ncoeff,ncoeff2,ncoeff3,ncoeff_max,nu
1448  integer :: nsym,shift_atm1(3)
1449  integer :: shift_atm2(3)
1450  real(dp):: dist_orig,dist_sym,tolsym8
1451  logical :: found,check_pbc,possible
1452 !arrays
1453  integer :: isym_rec(3,3),isym_rel(3,3),sc_size_in(3)
1454  integer :: transl(3),min_range(3),max_range(3)
1455  integer,allocatable :: blkval(:,:,:,:,:),list(:),list_symcoeff_tmp(:,:,:),list_symcoeff_tmp2(:,:,:)
1456  integer,allocatable :: list_symstr_tmp(:,:,:),indsym(:,:,:) ,symrec(:,:,:),symrel(:,:,:),list_symcoeff_tmp3(:,:,:)
1457  integer,allocatable :: index_irred(:)
1458  real(dp),allocatable :: tnons(:,:)
1459  real(dp),allocatable :: wkdist(:),distance(:,:,:)
1460  real(dp) :: difmin(3)
1461  real(dp) :: rprimd(3,3)
1462  real(dp) :: tratom(3)
1463  character(len=500) :: message
1464 
1465 
1466 
1467 !Initialisation of variables
1468  irpt_sym = 0
1469  nsym   = crystal%nsym
1470  rprimd = crystal%rprimd
1471  !ABI_MALLOC(xcart,(3,natom))
1472  !ABI_MALLOC(xred,(3,natom))
1473  !xcart(:,:) = crystal%xcart(:,:)
1474  !xred(:,:)  = crystal%xred(:,:)
1475  ncoeff_max = nrpt*natom*natom*3*3
1476 
1477 !Found the ref cell
1478  irpt_ref = 1
1479  do irpt=1,nrpt
1480    if(all(cell(:,irpt)==0))then
1481      irpt_ref = irpt
1482 !     exit
1483    end if
1484  end do
1485 
1486  !Set the size of the interaction
1487  check_pbc = .FALSE.
1488  sc_size_in = 0
1489  min_range = 0; max_range = 0
1490  if(present(sc_size))then
1491    sc_size_in = sc_size
1492    do mu=1,3
1493      call findBound_supercell(min_range(mu),max_range(mu),sc_size_in(mu))
1494    end do
1495  end if
1496 
1497  !Check which atom to fit, if not present do all atoms
1498  if(present(fit_iatom))then
1499     fit_iatom_in = fit_iatom
1500  else
1501     fit_iatom_in = -1
1502  endif
1503  if(fit_iatom_in > natom)then
1504       write(message, '(3a)' )&
1505 &       'fit_iatom cannot be greater than the number of atoms on the reference unit cell',ch10,&
1506 &       'Action: Change input'
1507    ABI_ERROR(message)
1508  end if
1509 
1510 !Obtain a list of rotated atom labels:
1511  ABI_MALLOC(indsym,(4,nsym,natom))
1512  ABI_MALLOC(symrec,(3,3,nsym))
1513  ABI_MALLOC(symrel,(3,3,nsym))
1514  ABI_MALLOC(tnons,(3,nsym))
1515  symrec = crystal%symrec
1516  symrel = crystal%symrel
1517  tnons  = crystal%tnons
1518 
1519  tolsym8=tol13
1520  call symatm(indsym,natom,nsym,symrec,tnons,&
1521 &            tolsym8,crystal%typat,crystal%xred)
1522  ABI_MALLOC(blkval,(3,natom,3,natom,nrpt))
1523  ABI_MALLOC(list,(natom*nrpt))
1524  ABI_MALLOC(list_symcoeff_tmp,(5,ncoeff_max,nsym))
1525  ABI_MALLOC(wkdist,(natom*nrpt))
1526 
1527 !1-Fill strain list
1528  ABI_MALLOC(list_symstr_tmp,(6,nsym,2))
1529  list_symstr_tmp = 1
1530  do ia=1,6
1531    if(list_symstr_tmp(ia,1,1)==0)cycle
1532 !  Transform the voigt notation
1533    if(ia<=3)then
1534      mu=ia;nu=ia
1535    else
1536      select case(ia)
1537      case(4)
1538        mu=2;nu=3
1539      case(5)
1540        mu=1;nu=3
1541      case(6)
1542        mu=1;nu=2
1543      end select
1544    end if
1545    do isym=1,nsym
1546 !    Get the symmetry matrix
1547      isym_rel(:,:) = crystal%symrel(:,:,isym)
1548      do idisy1=1,3
1549        do idisy2=1,3
1550          if((isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)) then
1551 !          Transform to the voig notation
1552            if(idisy1==idisy2)then
1553              list_symstr_tmp(ia,isym,1) = idisy1
1554              list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1)
1555            else
1556              if(idisy1==1.or.idisy2==1)then
1557                if(idisy1==2.or.idisy2==2)then
1558                  list_symstr_tmp(ia,isym,1) = 6
1559                end if
1560                if(idisy1==3.or.idisy2==3)then
1561                  list_symstr_tmp(ia,isym,1) = 5
1562                end if
1563              else
1564                list_symstr_tmp(ia,isym,1) = 4
1565              end if
1566            end if
1567            list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1) * isym_rel(nu,idisy2)
1568          end if
1569        end do
1570      end do
1571 !    Remove the symetric
1572 !     if(list_symstr_tmp(ia,isym,1) > ia) then
1573 !       list_symstr_tmp(list_symstr_tmp(ia,isym,1),:,1) = 0
1574 !     end if
1575    end do
1576  end do
1577 
1578 !Count the number of strain and transfert into the final array
1579   nstr_sym = 0
1580   do ia=1,6
1581     if(list_symstr_tmp(ia,1,1)/=0) nstr_sym = nstr_sym + 1
1582   end do
1583 
1584  if(allocated(list_symstr))then
1585    ABI_FREE(list_symstr)
1586  end if
1587  ABI_MALLOC(list_symstr,(nstr_sym,nsym,2))
1588 
1589  icoeff_tmp = 1
1590  do ia=1,6
1591    if(list_symstr_tmp(ia,1,1)/=0) then
1592      list_symstr(icoeff_tmp,:,:) = list_symstr_tmp(ia,:,:)
1593      icoeff_tmp = icoeff_tmp + 1
1594    end if
1595  end do
1596 !END STRAIN
1597 
1598 !Compute the distance between each atoms. Indeed the dist array contains the difference of
1599 !cartesian coordinate for each direction
1600 ! dist: vector between atom a, (0,0,0) and atom b (rpt)
1601 
1602  ABI_MALLOC(distance,(natom,natom,nrpt))
1603  ! Fortran 2008:  distance = norm2(dist, dim=1)
1604  do ia=1,natom
1605    do ib=1,natom
1606      do irpt=1,nrpt
1607        distance(ia,ib,irpt) = ((dist(1,ia,ib,irpt))**2+(dist(2,ia,ib,irpt))**2+&
1608 &                              (dist(3,ia,ib,irpt))**2)**0.5
1609      end do
1610    end do
1611  end do
1612 
1613 
1614 !Set to one blkval, all the coeff have to be compute
1615  blkval = 1
1616  icoeff = 1
1617  icoeff_tot = 1
1618  list_symcoeff_tmp = 0
1619 
1620 !2-Fill atom list
1621 !Big loop over generic atom
1622  do ia=1,natom
1623    wkdist(:)=reshape(distance(ia,:,:),(/natom*nrpt/))
1624    do ii=1,natom*nrpt
1625      list(ii)=ii
1626    end do
1627    ! FIXME: hexu: I think this should be improved.
1628    ! It seems that it depends on the specific order of the
1629    ! cell, and the implementation of sort algorithm.
1630    ! The sort should preserve the order if two distances are the same.
1631    ! In the future, the order of the cell might need to be unified
1632    ! Also it seems to depend on the cubic cell.
1633    call sort_dp(natom*nrpt,wkdist,list,tol8)
1634    do ii=1,natom*nrpt
1635 !    Get the irpt and ib
1636      irpt=(list(ii)-1)/natom+1
1637      ib=list(ii)-natom*(irpt-1)
1638      possible = .true.
1639 !Old way with the cut off
1640      if(cutoff < distance(ia,ib,irpt))then
1641        possible = .false.
1642      else
1643        ! in each direction jj, the d component d(jj) should be smaller than range_ifc(jj)
1644        do jj=1,3
1645 !          if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)  > tol10.or.&
1646 ! &            abs(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj))  < tol10)then
1647           if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)  > tol10)then
1648          possible = .false.
1649          end if
1650        end do
1651     endif
1652 
1653 !    If this distance is superior to the cutoff, we don't compute that term
1654      if(.not.possible)then
1655        blkval(:,ia,:,ib,irpt)= 0
1656        ! remove duplication: if in reference cell, only consider ab, not ba.
1657        if(irpt==irpt_ref)blkval(:,ib,:,ia,irpt)= 0
1658 !        Stop the loop
1659        cycle
1660      end if
1661 
1662 !    If this coefficient is not possible, we cycle...
1663      if (all(blkval(:,ia,:,ib,irpt)==0)) cycle
1664 
1665 !    Save the distance between the two atoms for futur checks
1666      dist_orig = (dist(1,ia,ib,irpt)**2+dist(2,ia,ib,irpt)**2+dist(3,ia,ib,irpt)**2)**0.5
1667 
1668      do mu=1,3
1669        do nu=1,3
1670 !      Check if : - The coefficient is not yet compute
1671 !                 - The directions are the same
1672 !                 - The atoms are not equivalent
1673          if (mu/=nu) then
1674            blkval(mu,ia,nu,ib,irpt)=0
1675            blkval(nu,ia,mu,ib,irpt)=0
1676            cycle
1677          end if
1678 !        Pass if the atoms are identical and in the ref cell
1679          if(irpt==irpt_ref.and.ia==ib) then
1680            blkval(mu,ia,nu,ib,irpt)=0
1681            blkval(nu,ib,mu,ia,irpt)=0
1682            cycle
1683          end if
1684          if(blkval(mu,ia,nu,ib,irpt)==1)then
1685 !          Loop over symmetries
1686            do isym=1,nsym
1687 !            Get the symmetry matrix for this sym
1688              isym_rec(:,:)  = crystal%symrec(:,:,isym)
1689              isym_rel(:,:) = crystal%symrel(:,:,isym)
1690 !            Get the corresponding atom and shift with the symetries
1691 !            For atom 1
1692              ipesy1 = indsym(4,isym,ia)
1693              shift_atm1 = indsym(1:3,isym,ia)
1694 !            And atom 2
1695              do jj=1,3 ! Apply transformation to original coordinates.
1696                tratom(jj) = dble(isym_rec(1,jj))*(crystal%xred(1,ib)+cell(1,irpt)-tnons(1,isym))&
1697                  &         +dble(isym_rec(2,jj))*(crystal%xred(2,ib)+cell(2,irpt)-tnons(2,isym))&
1698                  &         +dble(isym_rec(3,jj))*(crystal%xred(3,ib)+cell(3,irpt)-tnons(3,isym))
1699 
1700              end do
1701 
1702 !            Find symmetrically equivalent atom
1703              call symchk(difmin,ipesy2,natom,tratom,transl,crystal%typat(ib),&
1704 &                        crystal%typat,crystal%xred(:,:))
1705 
1706 !            Put information into array indsym: translations and label
1707              shift_atm2(:)= transl(:) - shift_atm1(:)
1708              found = .false.
1709              do irpt3=1,nrpt
1710                if(cell(1,irpt3)==shift_atm2(1).and.&
1711 &                 cell(2,irpt3)==shift_atm2(2).and.&
1712 &                 cell(3,irpt3)==shift_atm2(3))then
1713                  found = .true.
1714                  irpt_sym = irpt3
1715                end if
1716              end do
1717 
1718 !            Check the distance
1719              dist_sym = (dist(1,ipesy1,ipesy2,irpt_sym)**2+&
1720 &                        dist(2,ipesy1,ipesy2,irpt_sym)**2+&
1721 &                        dist(3,ipesy1,ipesy2,irpt_sym)**2)**0.5
1722              if(abs(dist_orig - dist_sym) > tol10)then
1723                write(message, '(a,i0,2a,I0,a,es15.8,2a,es15.8,2a)' )&
1724 &                'The distance between the atoms for the coefficient number ',icoeff,ch10,&
1725 &                'with the symmetry ',isym,' is ',dist_sym,ch10,'but the original distance is',&
1726 &                   dist_orig,ch10,&
1727 &                'Action: Contact abinit group'
1728                ABI_BUG(message)
1729              end if
1730 !            Now that a symmetric perturbation has been obtained,
1731 !            including the expression of the symmetry matrix, see
1732 !            if the symmetric perturbations are available
1733              do idisy1=1,3
1734                do idisy2=1,3
1735                  if (idisy1/=idisy2) then
1736 !                  Remove this term (is not computed)
1737 !                  Also remove opposite term... (Srx-Tix) = (Tix-Srx)
1738                    blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1739                    blkval(idisy2,ipesy1,idisy1,ipesy2,irpt_sym) = 0
1740                    cycle
1741                  else
1742                    if(isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)then
1743                      if(.not.found.or.(irpt_sym==irpt_ref.and.ipesy1==ipesy2)) then
1744 !                      Remove this term (is not computed) Sr-Sr or not include in the cell
1745 !                      Also remove oposite term... (Srx-Tix) = (Ti-Srx)
1746                        blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1747                        blkval(idisy2,ipesy2,idisy1,ipesy1,irpt_sym) = 0
1748                        cycle
1749                      else
1750 !                      Fill the list with the coeff and symmetric (need all symmetrics)
1751                        list_symcoeff_tmp(1:4,icoeff,isym)=(/idisy1,ipesy1,ipesy2,irpt_sym/)
1752 !                      Check the sign
1753                        if(isym_rel(mu,idisy1)/=isym_rel(nu,idisy2))then
1754                          write(message, '(a,i0,a,I0,4a)' )&
1755 &                        'The sign of coefficient number ',icoeff,' with the symmetry ',isym,ch10,&
1756 &                        'can not be found... Something is going wrong',ch10,&
1757 &                        'Action: Contact abinit group'
1758                          ABI_BUG(message)
1759                        end if
1760                        list_symcoeff_tmp(5,icoeff,isym)= isym_rel(nu,idisy2)
1761                      end if
1762                    end if
1763                  end if
1764                end do
1765              end do
1766            end do ! end loop sym
1767            icoeff = icoeff + 1
1768          end if
1769 !        This coeff is now computed
1770          blkval(mu,ia,nu,ib,irpt)= 0
1771        end do ! end loop nu
1772      end do ! end loop mu
1773    end do ! end loop ii: iatom and irpt.
1774  end do ! end loop ia
1775 
1776 !Reset the output
1777  if(allocated(list_symcoeff))then
1778    ABI_FREE(list_symcoeff)
1779  end if
1780 
1781  ABI_FREE(distance)
1782 
1783 !Transfert the final array with all the coefficients
1784 !With this array, we can access to all the terms presents
1785 !ncoeff1 + symetrics
1786 !first dimension is 4 (mu,ia,ib,irpt)
1787 !irpt is the index of the cell of the atom ib in the cell array
1788 !example cell(:,irpt=12) can be (-1 0 -2). The cell of ia is
1789 !always 0 0 0
1790 !Transfert the final array for the list of irreductible coeff and symetries
1791 !With this array, we can access to the irretuctible coefficients (ncoeff1) and
1792 !all the symetrics of these coefficients (nsym)
1793 !first dimension is 5 (mu,ia,ib,irpt,icoeff)
1794 !icoeff is the position of this coefficients in the list_fullcoeff array
1795 
1796 !1/ step remove the zero coeff in this array
1797  ncoeff = 0
1798  do icoeff = 1,ncoeff_max
1799      if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
1800      ncoeff = ncoeff + 1
1801    end if
1802  end do
1803 
1804  ABI_MALLOC(list_symcoeff_tmp2,(6,ncoeff,nsym))
1805  list_symcoeff_tmp2 = 0
1806  icoeff = 0
1807  do icoeff_tmp = 1,ncoeff_max
1808    if(.not.(all(list_symcoeff_tmp(:,icoeff_tmp,1)==0)))then
1809      icoeff = icoeff + 1
1810      list_symcoeff_tmp2(1:5,icoeff,:) = list_symcoeff_tmp(1:5,icoeff_tmp,:)
1811      end if
1812  end do
1813 
1814 
1815 !2/ set the dimension six of list_symcoeff_tmp2(6,icoeffs,1)
1816 !   and check is a symetric coeff is not coresspondig to an other
1817 !   one, in this case we set this coeff to 0
1818 ! ncoeff2 = zero
1819  do icoeff = 1,ncoeff
1820 !  found the index of each coeff in list_fullcoeff
1821    do isym = 1,nsym
1822      icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,1),&
1823 &                               list_symcoeff_tmp2(2,icoeff,isym),&
1824 &                               list_symcoeff_tmp2(3,icoeff,isym),&
1825 &                               list_symcoeff_tmp2(4,icoeff,isym),&
1826 &                               list_symcoeff_tmp2(1,icoeff,isym),&
1827 &                               ncoeff)
1828      list_symcoeff_tmp2(6,icoeff,isym) = icoeff2
1829    end do
1830  end do
1831 
1832 !2.5/do checks
1833  do icoeff = 1,ncoeff
1834    do isym = 1,nsym
1835      if(list_symcoeff_tmp2(6,icoeff,isym)==0)then
1836        write(message, '(a,i0,a,I0,4a)' )&
1837 &           'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1838 &           'has no equivalent',ch10,&
1839 &           'Action: Contact abinit group'
1840        ABI_BUG(message)
1841      else
1842        if(icoeff /= list_symcoeff_tmp2(6,icoeff,isym))then
1843          if(list_symcoeff_tmp2(1,icoeff,isym)/=&
1844 &           list_symcoeff_tmp2(1,list_symcoeff_tmp2(6,icoeff,isym),1))then
1845            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1846 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1847 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1848 &          'because the direction is different:',ch10,&
1849 &          'Action: Contact abinit group'
1850            ABI_BUG(message)
1851          end if
1852          if(list_symcoeff_tmp2(4,icoeff,isym)/=&
1853 &           list_symcoeff_tmp2(4,list_symcoeff_tmp2(6,icoeff,isym),1))then
1854            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1855 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1856 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1857 &          'because the cell is different',ch10,&
1858 &          'Action: Contact abinit group'
1859            ABI_BUG(message)
1860          end if
1861          if((list_symcoeff_tmp2(2,icoeff,isym)/=&
1862 &            list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1863 &            list_symcoeff_tmp2(3,icoeff,isym)/=&
1864 &            list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1)).and.&
1865 &           (list_symcoeff_tmp2(2,icoeff,isym)/=&
1866 &            list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1867 &            list_symcoeff_tmp2(3,icoeff,isym)/=&
1868 &            list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1)))then
1869            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1870 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1871 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1872 &          'because the atoms different',ch10,&
1873 &          'Action: Contact abinit group'
1874            ABI_BUG(message)
1875          end if
1876        end if
1877      end if
1878    end do
1879  end do
1880 
1881 
1882 !Check if the atom 2 is not (with the PBC) is in the same cell than
1883 !the atom 1. For example if you want to fit 2x2x2 cell the interation
1884 !Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process...
1885  do icoeff = 1,ncoeff
1886    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1887      do mu=1,3
1888        ! out of max_range(mu)
1889        if( -max_range(mu) > cell(mu,list_symcoeff_tmp2(4,icoeff,1)) .or. &
1890 &           cell(mu,list_symcoeff_tmp2(4,icoeff,1)) > max_range(mu))then
1891 
1892          ! FIXME: hexu: cell(mu, irpt). But why -1?
1893          ! The if above looks sufficient.
1894          !if(cell(mu,list_symcoeff_tmp2(4,icoeff,1)) < -1)then
1895          list_symcoeff_tmp2(:,icoeff,:)=0
1896          exit
1897          !end if
1898        end if
1899      end do
1900    end if
1901  !MS only keep terms with ia == fit_iatom or ib == fit_iatom specified in input
1902    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. fit_iatom_in > 0)then
1903      if (list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then !& LB
1904 !&        list_symcoeff_tmp2(3,icoeff,1) /= fit_iatom_in)then  LB
1905          list_symcoeff_tmp2(:,icoeff,1) = 0
1906 !     else if(list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then ! &
1907 !&            .and. list_symcoeff_tmp2(3,icoeff,1) == fit_iatom_in &  !LB
1908 !&            .and. any(cell(:,list_symcoeff_tmp2(4,icoeff,1)) /= 0))then !LB
1909 !         list_symcoeff_tmp2(:,icoeff,1) = 0
1910      endif
1911    endif
1912 !For Debugging keep only terms A_x-A_x[100] etc. comment if above
1913 !  if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1914 !    if (list_symcoeff_tmp2(2,icoeff,1) /= 1 .or. list_symcoeff_tmp2(3,icoeff,1) /= 1)then
1915 !        list_symcoeff_tmp2(:,icoeff,1) = 0
1916 !    endif  !MS only keep terms with ia == 1
1917 !  endif
1918 end do
1919 
1920 !write(std_out,*) "DEBUG: min_range(:): ", min_range(:)
1921 !write(std_out,*) "DEBUG: max_range(:): ", max_range(:)
1922 
1923 !3/ Remove useless terms like opposites
1924  do icoeff = 1,ncoeff
1925    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. & ! valid term.
1926 &      list_symcoeff_tmp2(2,icoeff,1) /= list_symcoeff_tmp2(3,icoeff,1))then ! on-site term. ! FIXME: hexu:I don't understand if iatom=jatom but rpt/=0, it is not onsite.
1927    do isym = 1,nsym
1928    !icoeff2 = list_symcoeff_tmp2(6,icoeff,isym)
1929    !if (icoeff2> icoeff)then
1930    !  list_symcoeff_tmp2(:,icoeff2,1) = 0
1931    !end if
1932      do jsym=1,nsym
1933        ! FIXME: Should the rpt be -rpt??
1934 !       icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),&
1935 !&                                 list_symcoeff_tmp2(3,icoeff,isym),&
1936 !&                                 list_symcoeff_tmp2(2,icoeff,isym),&
1937 !&                                 list_symcoeff_tmp2(4,icoeff,isym),& !rpt
1938 !&                                 list_symcoeff_tmp2(1,icoeff,isym),&
1939 !&                                 ncoeff)
1940        icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),&
1941 &                                 list_symcoeff_tmp2(3,icoeff,isym),&
1942 &                                 list_symcoeff_tmp2(2,icoeff,isym),&
1943 &                                 list_symcoeff_tmp2(4,icoeff,isym),& !rpt
1944 &                                 find_opposite_irpt(cell, list_symcoeff_tmp2(1,icoeff,isym)),&
1945 &                                 ncoeff)
1946 
1947        if (icoeff2> icoeff)then
1948          list_symcoeff_tmp2(:,icoeff2,1) = 0
1949        end if
1950      end do
1951    end do
1952    end if
1953  end do
1954 
1955 !4/ Recount the number of coeff after step 3
1956  ncoeff2 = 0
1957  isym = 0
1958  do icoeff = 1,ncoeff
1959    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1960      ncoeff2 = ncoeff2 + 1
1961    end if
1962  end do
1963 
1964 !write(std_out,*) "DEBUG ncoeff after 2.5.3: ", ncoeff
1965 !write(std_out,*) "DEBUG ncoeff2 after 2.5.3: ", ncoeff2
1966 
1967  ABI_FREE(list_symcoeff_tmp)
1968  ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff,nsym))
1969  list_symcoeff_tmp = list_symcoeff_tmp2
1970 
1971 !4.1 Count irreducible terms
1972  do icoeff = 1,ncoeff
1973    do isym = 1,nsym
1974      icoeff2 = list_symcoeff_tmp2(6,icoeff,isym)
1975      if (icoeff2> icoeff)then
1976        list_symcoeff_tmp(:,icoeff2,1) = 0
1977      end if
1978    end do
1979  end do
1980 
1981 ncoeff3 = 0
1982  do icoeff = 1,ncoeff
1983    if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
1984      ncoeff3 = ncoeff3 + 1
1985    end if
1986  end do
1987 !write(std_out,*) "DEBUG ncoeff3 after 2.5.3: ", ncoeff3
1988 
1989 !4.2 Put irreducible terms in front of list_symcoeff_tmp3
1990 !    Store index of irreducible terms
1991  ABI_MALLOC(list_symcoeff_tmp3,(6,ncoeff2,nsym))
1992  ABI_MALLOC(index_irred,(ncoeff3))
1993 icoeff_tmp = 0
1994 do icoeff = 1,ncoeff
1995    if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
1996      icoeff_tmp = icoeff_tmp+1
1997      index_irred(icoeff_tmp) = icoeff
1998      list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp(:,icoeff,:)
1999    endif
2000 enddo
2001 
2002 !4.3 Put symmetric equivalents behind in list_symcoeff_tmp3
2003 !    Attention icoeff_tmps keeps it's value of loop before
2004 !    TODO for check should be equal to ncoeff2
2005 do icoeff = 1,ncoeff
2006    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0))&
2007 &     .and..not. any(index_irred == icoeff))then
2008      icoeff_tmp = icoeff_tmp+1
2009      list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp2(:,icoeff,:)
2010    endif
2011 enddo
2012 
2013 !4.4 A little copy round
2014 ABI_FREE(list_symcoeff_tmp)
2015  ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff2,nsym))
2016 list_symcoeff_tmp = list_symcoeff_tmp3
2017 
2018 !5/ Final transfert
2019  ABI_MALLOC(list_symcoeff,(6,ncoeff2,nsym))
2020  list_symcoeff = 0
2021  icoeff = 0
2022  do icoeff = 1,ncoeff2
2023    list_symcoeff(1:6,icoeff,:) = list_symcoeff_tmp(1:6,icoeff,:)
2024    do isym=1,nsym
2025    end do
2026  end do
2027 
2028 !6/ reset the dimension six of list_symcoeff_tmp2(6,icoeffs,1)
2029 !   and check is a symetric coeff is not coresspondig to an other
2030 !   one, in this case we set this coeff to 0
2031  do icoeff = 1,ncoeff2
2032 !  found the index of each coeff in list_fullcoeff
2033    do isym = 1,nsym
2034      icoeff2 = getCoeffFromList(list_symcoeff(:,:,1),&
2035 &                               list_symcoeff(2,icoeff,isym),&
2036 &                               list_symcoeff(3,icoeff,isym),&
2037 &                               list_symcoeff(4,icoeff,isym),&
2038 &                               list_symcoeff(1,icoeff,isym),&
2039 &                               ncoeff2)
2040      list_symcoeff(6,icoeff,isym) = icoeff2
2041 !     list_symcoeff(6,icoeff2,isym) = icoeff
2042    end do
2043  end do
2044 
2045 !Set the max number of coeff inside list_symcoeff
2046  ncoeff_sym = ncoeff3
2047 
2048 
2049 !Deallocation
2050  ABI_FREE(blkval)
2051  ABI_FREE(list)
2052  ABI_FREE(list_symcoeff_tmp)
2053  ABI_FREE(list_symcoeff_tmp2)
2054  ABI_FREE(list_symcoeff_tmp3)
2055  ABI_FREE(list_symstr_tmp)
2056  ABI_FREE(index_irred)
2057  ABI_FREE(indsym)
2058  ABI_FREE(symrec)
2059  ABI_FREE(symrel)
2060  ABI_FREE(tnons)
2061  !ABI_FREE(xcart)
2062  !ABI_FREE(xred )
2063  ABI_FREE(wkdist)
2064 end subroutine polynomial_coeff_getList

m_polynomial_coeff/polynomial_coeff_getName [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_getName

FUNCTION

 get the name of a polynomial coefficient

INPUTS

 natom = number of atoms
 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
 symbols(natom)  =  array with the atomic symbol:["Sr","Ru","O1","O2","O3"]
 recompute = (optional) flag to set if the name has to be recomputed
 iterm = (optional) number of the term used for the name

OUTPUT

 name = name xof the coefficients

SOURCE

473 subroutine polynomial_coeff_getName(name,polynomial_coeff,symbols,recompute,iterm)
474 
475  implicit none
476 
477 !Arguments ------------------------------------
478 !scalars
479  integer,optional,intent(in) :: iterm
480 !arrays
481  character(len=5),intent(in) :: symbols(:)
482  character(len=200),intent(out):: name
483  type(polynomial_coeff_type),optional, intent(in) :: polynomial_coeff
484  logical,optional,intent(in) :: recompute
485 !Local variables-------------------------------
486 !scalar
487  integer :: ii,idisp,iterm_in
488  logical :: need_recompute
489 !arrays
490  integer :: cell_atm1(3),cell_atm2(3)
491  character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
492  character(len=1) :: dir
493  character(len=2) :: power_disp,power_dispchar
494  character(len=20) :: atm1,atm2
495  character(len=100):: atm1_tmp,atm2_tmp
496  character(len=200):: text
497  character(len=500):: message
498  ! *************************************************************************
499 
500  !Reset output
501  name=""
502  iterm_in = 1
503 
504  !Set the optional arguments
505  need_recompute = .FALSE.
506  if(present(recompute)) need_recompute = recompute
507  if(present(iterm)) then
508    iterm_in = iterm
509  else
510    if(need_recompute)then
511      iterm_in = -1
512      do ii=1,polynomial_coeff%nterm
513        !      Find the index of the ref
514        if(iterm_in==-1) then !Need to find the reference term
515          do idisp=1,polynomial_coeff%terms(ii)%ndisp
516            if(polynomial_coeff%terms(ii)%direction(idisp) > 0) then
517              iterm_in = ii
518              if(any(polynomial_coeff%terms(ii)%cell(:,1,idisp) /= 0).or.&
519                &               any(polynomial_coeff%terms(ii)%cell(:,2,idisp) /= 0)) then
520                iterm_in = -1
521                exit
522              end if
523            end if
524          end do!end do disp
525        else
526          exit
527        end if
528      end do!end do term
529      !    If not find, we set to the first element
530      if(iterm_in==-1) iterm_in = 1
531    else
532      iterm_in = 1
533    end if
534  end if
535  !Do check
536  if(iterm_in > polynomial_coeff%nterm.or.iterm_in < 0) then
537    write(message, '(5a)')&
538      &      ' The number of the requested term for the generation of',ch10,&
539      &      'the name of the coefficient is not possible.',ch10,&
540      &      'Action: Contact Abinit group.'
541    ABI_BUG(message)
542  end if
543 
544  if(polynomial_coeff%name /= "".and..not.need_recompute)then
545    name = polynomial_coeff%name
546  else
547 !  Nedd to recompute
548    do idisp=1,polynomial_coeff%terms(iterm_in)%ndisp
549      text = ""
550      !Fill variables for this displacement
551      write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_disp(idisp)
552      power_disp=trim(power_dispchar)
553 
554      atm1=symbols(polynomial_coeff%terms(iterm_in)%atindx(1,idisp))
555      atm2=symbols(polynomial_coeff%terms(iterm_in)%atindx(2,idisp))
556      dir=mutodir(polynomial_coeff%terms(iterm_in)%direction(idisp))
557      cell_atm1=polynomial_coeff%terms(iterm_in)%cell(:,1,idisp)
558      cell_atm2=polynomial_coeff%terms(iterm_in)%cell(:,2,idisp)
559 !    Construct ATM1
560      if (any(cell_atm1(:) /= 0) )then
561        write(atm1_tmp,'(4a,I0,a,I0,a,I0,a)')  trim(atm1),"_",dir,"[",cell_atm1(1)," ",&
562 &                                               cell_atm1(2)," ",cell_atm1(3),"]"
563        else
564        atm1_tmp = trim(atm1)//"_"//dir
565      end if
566 !      Construct ATM2
567      if(any(cell_atm2(:) /= 0))then
568        write(atm2_tmp,'(4a,I0,a,I0,a,I0,a)')  trim(atm2),"_",dir,"[",cell_atm2(1)," ",&
569  &                                              cell_atm2(2)," ",cell_atm2(3),"]"
570      else
571        atm2_tmp = trim(atm2)//"_"//dir
572      end if
573      text="("//trim(atm1_tmp)//"-"//trim(atm2_tmp)//")^"//power_disp
574      name = trim(name)//trim(text)
575    end do
576    !Strain case
577    do idisp=1,polynomial_coeff%terms(iterm_in)%nstrain
578      write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_strain(idisp)
579      power_disp=trim(power_dispchar)
580      dir=mutodir(3+polynomial_coeff%terms(iterm_in)%strain(idisp))
581      text="("//"eta_"//trim(dir)//")^"//power_disp
582      name = trim(name)//trim(text)
583    end do
584  end if
585 
586 end subroutine polynomial_coeff_getName

m_polynomial_coeff/polynomial_coeff_getNorder [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_getNorder

FUNCTION

 Compute and store into the datatype coefficients, all the possible
 coefficients for given orders

INPUTS

 cutoff = cut-off for the inter atomic forces constants
 crystal<type(crystal_t)> = datatype with all the information for the crystal
 power_disps(2) = array with the minimal and maximal power_disp to be computed
 max_power_strain = maximum order of the strain of the strain phonon coupling
 option = 0 compute all terms
          1 still in development
 sc_size(3) = size of the supercell used for the fit.
               For example if you want to fit 2x2x2 cell the interation
               Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process
               If check_pbc is true we remove these kind of terms
 comm = MPI communicator
 anharmstr = logical, optional : TRUE, the anharmonic strain is computed (\eta)^power_disp ...
                                   FALSE, (default) the anharmonic strain are not computed
 spcoupling= logical, optional : TRUE(default) the anharmonic strain-phonon coupling is computed
 only_odd_power = logical, optional : if TRUE return only odd power
 only_even_power= logical, optional : if TRUe return only even power
 distributed = logical, optional : True, the coefficients will be distributed on the CPU
 verbose  = optional, flag for the verbose mode

OUTPUT

 polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff) = array of datatype with the polynomial_coeff
 ncoeff = number of coefficients for this CPU if distributed == true, all otherwise
 ncoeff_tot = total number of coefficient over the CPU

SOURCE

2104 subroutine polynomial_coeff_getNorder(coefficients,crystal,cutoff,ncoeff,ncoeff_tot,power_disps,&
2105 &                                     max_power_strain,option,sc_size,comm,anharmstr,spcoupling,&
2106 &                                     distributed,only_odd_power,only_even_power,fit_iatom,&
2107 &                                     compute_symmetric,dispterms,verbose)
2108 
2109  implicit none
2110 
2111 !Arguments ------------------------------------
2112 !scalars
2113  integer,intent(in) :: max_power_strain,option,comm
2114  integer,intent(out):: ncoeff,ncoeff_tot
2115  integer,optional,intent(in) :: fit_iatom
2116  real(dp),intent(in):: cutoff
2117  logical,optional,intent(in) :: anharmstr,spcoupling,distributed,verbose,dispterms
2118  logical,optional,intent(in) :: only_odd_power,only_even_power,compute_symmetric
2119 !arrays
2120  integer,intent(in) :: power_disps(2),sc_size(3)
2121  type(crystal_t), intent(inout) :: crystal
2122  type(polynomial_coeff_type),allocatable,intent(inout) :: coefficients(:)
2123 !Local variables-------------------------------
2124 !scalar
2125  integer :: icoeff,icoeff2,icoeff3,ierr,ii,iterm
2126  integer :: i
2127  integer :: master,my_rank,my_ncoeff,my_newncoeff,natom,ncombination,ncoeff_max,ncoeff_sym
2128  integer :: ncoeff_symsym,nirred_comb,iirred_comb,ndisp,nstrain,fit_iatom_in
2129  integer :: ncoeff_alone,ndisp_max,nproc,nrpt,nsym,nterm,nstr_sym,my_size
2130  integer :: my_icoeff,rank_to_send,rank_to_receive,rank_to_send_save
2131  integer :: ncombi_alone,my_ncombi_simple,my_ncombi_start,my_ncombi_end,my_ncombi,my_nirred
2132  logical :: iam_master,need_anharmstr,need_spcoupling,need_distributed,need_verbose
2133  logical :: need_only_odd_power,need_only_even_power,compute_sym,irreducible,need_compute_symmetric
2134  logical :: need_dispterms
2135 !arrays
2136  integer :: shape_listsymcoeff(3),shape_listsymstr(3)
2137  integer,allocatable :: buffsize(:),buffdispl(:) !,dummylist(:) ,index_irred(:)
2138  integer,allocatable :: offsets(:)
2139  integer,allocatable :: cell(:,:),compatibleCoeffs(:,:)
2140  integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:),list_coeff(:),list_combination(:,:)
2141  integer,allocatable :: list_combination_tmp(:,:)
2142  integer,allocatable :: irank_ncombi(:),my_index_irredcomb(:)
2143  integer,allocatable :: my_coefflist(:),my_coeffindexes(:),my_newcoeffindexes(:),my_list_combination(:,:)
2144  type(int2d_array_type) :: my_array_combination
2145  integer,allocatable :: my_list_combination_tmp(:,:)
2146  real(dp) :: rprimd(3,3),range_ifc(3)
2147  real(dp),allocatable :: dist(:,:,:,:)
2148  character(len=5),allocatable :: symbols(:)
2149  character(len=200):: name
2150  character(len=500) :: message
2151  type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_tmp
2152  type(polynomial_term_type),dimension(:),allocatable :: terms
2153  character(len=fnlen) :: filename
2154 ! *************************************************************************
2155  !Hide filename for debugging
2156  ABI_UNUSED(filename)
2157 
2158  !MPI variables
2159  master = 0
2160  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2161  iam_master = (my_rank == master)
2162 
2163 !Free the output
2164  if(allocated(coefficients))then
2165    do ii =1,size(coefficients)
2166      call polynomial_coeff_free(coefficients(ii))
2167    end do
2168    ABI_FREE(coefficients)
2169  end if
2170 
2171 !Check
2172  if(option > power_disps(2))then
2173    write(message, '(3a)' )&
2174 &       'Option can not be superior to the maximum order ',ch10,&
2175 &       'Action: contact abinit group'
2176    ABI_ERROR(message)
2177  end if
2178 
2179 !Initialisation of variables
2180  need_anharmstr = .TRUE.
2181  if(present(anharmstr)) need_anharmstr = anharmstr
2182  need_spcoupling = .TRUE.
2183  if(present(spcoupling)) need_spcoupling = spcoupling
2184  need_distributed = .FALSE.
2185  if(present(distributed)) need_distributed  = distributed
2186  need_verbose = .TRUE.
2187  if(present(verbose)) need_verbose = verbose
2188  need_only_odd_power = .FALSE.
2189  if(present(only_odd_power)) need_only_odd_power = only_odd_power
2190  need_only_even_power = .FALSE.
2191  if(present(only_even_power)) need_only_even_power = only_even_power
2192  need_compute_symmetric = .TRUE.
2193  if(present(compute_symmetric)) need_compute_symmetric = compute_symmetric
2194  need_dispterms = .TRUE.
2195  if(present(dispterms)) need_dispterms = dispterms
2196 
2197  if(need_only_odd_power.and.need_only_even_power)then
2198       write(message, '(3a)' )&
2199 &       'need_only_odd_power and need_only_even_power are both true',ch10,&
2200 &       'Action: contact abinit group'
2201    ABI_ERROR(message)
2202  end if
2203 
2204  !Check which atom to fit, if not present do all atoms
2205  if(present(fit_iatom))then
2206    ! FIXME: the fit_iatom as the input is the index of the atom
2207    ! but in the fitting subroutine, it is the irreducible atom.
2208     fit_iatom_in = fit_iatom
2209  else
2210     fit_iatom_in = -1
2211  endif
2212 
2213   rprimd = crystal%rprimd
2214 
2215  !=
2216 call prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc , symbols)
2217 
2218 if(iam_master)then
2219  if(need_verbose)then
2220    write(message,'(1a)')' Generation of the list of all the possible pairs of atoms within cutoff'
2221    call wrtout(std_out,message,'COLL')
2222  end if
2223  call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
2224 &                              natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size,&
2225 &                              fit_iatom=fit_iatom_in)
2226 
2227 
2228  shape_listsymcoeff = shape(list_symcoeff)
2229  shape_listsymstr   = shape(list_symstr)
2230 endif!if iam master
2231 
2232 !Broadcast Results of getList
2233 call xmpi_bcast(shape_listsymcoeff, master, comm, ierr)
2234 call xmpi_bcast(shape_listsymstr, master, comm, ierr)
2235 call xmpi_bcast(nstr_sym, master, comm, ierr)
2236 call xmpi_bcast(ncoeff_sym, master, comm, ierr)
2237 if(.not. iam_master )then
2238   ABI_MALLOC(list_symcoeff,(shape_listsymcoeff(1),shape_listsymcoeff(2),shape_listsymcoeff(3)))
2239   ABI_MALLOC(list_symstr,(shape_listsymstr(1),shape_listsymstr(2),shape_listsymstr(3)))
2240 endif
2241 call xmpi_bcast(list_symcoeff, master, comm, ierr)
2242 call xmpi_bcast(list_symstr, master, comm, ierr)
2243 call xmpi_barrier(comm)
2244 
2245 ncoeff_symsym = size(list_symcoeff(1,:,1))
2246 
2247 !Compute the total number of coefficient
2248  ncoeff_tot = ncoeff_sym+nstr_sym
2249 
2250 !if(iam_master)then
2251 !Check the distanceance bewteen coefficients and store integer:
2252 ! 0: the mix between these coefficient is not possible
2253 ! 1: the mix between these coefficient is possible
2254    ABI_MALLOC(compatibleCoeffs,(ncoeff_symsym+nstr_sym,ncoeff_symsym+nstr_sym))
2255    compatibleCoeffs(:,:) = 1
2256 
2257    if(need_verbose)then
2258      write(message,'(1a)')' Check the compatible coefficients with respect to the cutoff'
2259      call wrtout(std_out,message,'COLL')
2260    end if
2261 
2262    do icoeff=1,ncoeff_symsym+nstr_sym
2263      do icoeff2=1,ncoeff_symsym+nstr_sym
2264 !      Select case:
2265 !      if both icoeff are displacement => check the distance
2266 !      if both icoeff are strain => check the flag
2267 !      Otherwise cycle (we keep the term)
2268 !       if(icoeff>ncoeff_sym.and.icoeff2<=ncoeff_sym)cycle
2269 !       if(icoeff<=ncoeff_sym.and.icoeff2>ncoeff_sym)cycle
2270        !Forbid anharmonic strain terms if not wanted
2271        if((icoeff>ncoeff_symsym .and. icoeff2>ncoeff_symsym).and.&
2272 &       .not.need_anharmstr) then
2273          compatibleCoeffs(icoeff,icoeff2) = 0
2274          compatibleCoeffs(icoeff2,icoeff) = 0
2275        end if
2276        !Forbid strain phonon terms if not wanted
2277        if((icoeff>ncoeff_symsym.or.icoeff2>ncoeff_symsym).and.&
2278 &       .not.need_spcoupling) then
2279          compatibleCoeffs(icoeff,icoeff2) = 0
2280          compatibleCoeffs(icoeff2,icoeff) = 0
2281        end if
2282 !!TEST_AM
2283 !           compatibleCoeffs(icoeff,icoeff2) = 0
2284 !           compatibleCoeffs(icoeff2,icoeff) = 0
2285 !!TEST_AM
2286 !         end if
2287 !       end if
2288 
2289        ! TODO: hexu: Check if this is only valid for orthogonal structures. As it takes the x component of the rprimd
2290        ! to compare.
2291        ! Here it checks: for two pairs of atom with icoeff (a-b)  and icoeff2 (c-d).
2292        !. (a) a-d along x. (b) a-d along y. (c) a-d along z are within the range.
2293        ! As a-b and c-d are within the cutoff, and a==c, there is no more need to check.
2294 
2295 
2296        if(icoeff<=ncoeff_symsym.and.icoeff2<=ncoeff_symsym)then !Check combination of irreducible bodies and their symmetric equivalent
2297          if(list_symcoeff(2,icoeff,1)/= list_symcoeff(2,icoeff2,1)) then
2298            ABI_BUG("The first components of the pairs in the coefficient are not equivalent.")
2299          end if
2300 
2301           if(abs(dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! rx(a, b)
2302 &            -dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! rx(a, d )
2303 &            >= (rprimd(1,1) + rprimd(1,2) + rprimd(1,3))*sc_size(1)  .or. &
2304              abs(dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) &  ! ry(a, b)
2305 &            -dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) &  ! ry(a, d)
2306 &            >= (rprimd(2,1) + rprimd(2,2) + rprimd(2,3))*sc_size(2)  .or. &
2307              abs(dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) &  ! rz(a, b)
2308 &            -dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) &  ! rz(a, d)
2309 &            >= (rprimd(3,1) + rprimd(3,2) + rprimd(3,3))*sc_size(3))then
2310              compatibleCoeffs(icoeff,icoeff2) = 0
2311              compatibleCoeffs(icoeff2,icoeff) = 0
2312           endif
2313        endif
2314      end do !end  icoeff
2315    end do !icoeff2
2316 
2317  ABI_FREE(dist)
2318 !  Compute all the combination of coefficient up to the given order  (get the number)
2319    if(need_verbose)then
2320      write(message,'(1a)')' Compute the number of possible combinations'
2321      call wrtout(std_out,message,'COLL')
2322    end if
2323 
2324    ABI_MALLOC(list_coeff,(0))
2325    ABI_MALLOC(list_combination,(0,0))
2326    icoeff  = 1
2327    icoeff2 = 0
2328    nirred_comb = 0
2329    iirred_comb = 0
2330    call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2331 &                   list_coeff,list_combination,icoeff,max_power_strain,natom,ncoeff_sym,&
2332 &                   ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,icoeff,nrpt,nsym,1,power_disps(1),power_disps(2),symbols,comm,&
2333 &                   nbody=option,compute=.false.,anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2334 &                   only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms)
2335    ABI_FREE(list_coeff)
2336    ABI_FREE(list_combination)
2337 
2338 !  Output how much we found
2339    if(need_verbose)then
2340      write(message,"(1a,I10)") "  -Number of irreducible pairs within cutoff:   ", ncoeff_sym
2341      call wrtout(std_out,message,'COLL')
2342      write(message,"(1a,I10)") "  -Number of combinations of irreducible pairs: ", nirred_comb
2343      call wrtout(std_out,message,'COLL')
2344      !write(message,"(1a,I10)") "  -Number of possible symmetric combinations:   ", ncombination
2345      !call wrtout(std_out,message,'COLL')
2346      write(message,'(a,a)') ch10,' Compute the combinations of irreducible pairs'
2347      call wrtout(std_out,message,'COLL')
2348    end if
2349    ABI_MALLOC(list_coeff,(0))
2350    ABI_MALLOC(list_combination_tmp,(power_disps(2),nirred_comb))
2351    icoeff  = 1
2352    icoeff2 = 0
2353    iirred_comb = 0
2354    list_combination_tmp = 0
2355 !  Compute all the combination of coefficient up to the given order
2356    call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2357 &                   list_coeff,list_combination_tmp,icoeff,max_power_strain,natom,&
2358 &                   ncoeff_sym,ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,ncombination,nrpt,nsym,1,power_disps(1),&
2359 &                   power_disps(2),symbols,comm,nbody=option,compute=.true.,&
2360 &                   anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2361 &                   only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms)
2362    ABI_FREE(list_coeff)
2363    nirred_comb = size(list_combination_tmp,2)
2364 
2365 ! If we want to compute equivalent symmetric combinations go here.
2366 if(need_compute_symmetric)then
2367   if(need_verbose)then
2368     write(message,'(1a)')' Distribute irreducible combinations over CPU'
2369     call wrtout(std_out,message,'COLL')
2370     write(message,'(1a)')' Compute symmetric combinations of combinations of irreducible pairs'
2371     call wrtout(std_out,message,'COLL')
2372     write(message,'(3a)')' ---> Try to match number of CPU to number combinations of irreducible pairs',ch10,&
2373     &                    '      for max. speedup'
2374     call wrtout(std_out,message,'COLL')
2375   endif
2376 
2377   ! COUNT IRREDUCIBLE COMBINATIONS FOR EACH PROCESSOR
2378   ncombi_alone = mod(nirred_comb,nproc)
2379   my_ncombi_simple = int(aint(real(nirred_comb,sp)/(nproc)))
2380   if(ncombi_alone == 0 .and. nirred_comb >= nproc)then !ncombi > nproc and no remainder
2381     my_ncombi_start = (my_ncombi_simple * my_rank) + 1
2382     my_ncombi_end   = my_ncombi_start + my_ncombi_simple - 1
2383   else if(nirred_comb < nproc)then  !ncombi smaller than nproc
2384     if(my_rank + 1 <= nirred_comb)then !myrank smaller than ncombi
2385       my_ncombi_start = my_rank + 1
2386       my_ncombi_end = my_ncombi_start
2387     else
2388       my_ncombi_start = nirred_comb + 1 !myrank bigger than ncombi
2389       my_ncombi_end = nirred_comb + 1
2390     endif
2391   else if(nirred_comb > nproc .and. ncombi_alone /= 0)then !ncombi > nproc and remainder
2392     if(my_rank >= (nproc-ncombi_alone)) then
2393       my_ncombi_start = (my_ncombi_simple * my_rank) + 1 + (my_rank - nproc + ncombi_alone)
2394       my_ncombi_end = my_ncombi_start + my_ncombi_simple
2395     else
2396       my_ncombi_start = (my_ncombi_simple * my_rank) + 1
2397       my_ncombi_end   = my_ncombi_start + my_ncombi_simple - 1
2398     endif
2399   end if
2400 
2401   if(my_ncombi_end <= nirred_comb)then
2402      my_nirred = my_ncombi_end+1-my_ncombi_start
2403   else
2404      my_nirred = 0
2405   endif
2406 
2407 
2408   !write(std_out,*) "my_rank", my_rank, "my_ncmobi_start", my_ncombi_start, "my_ncombi_end", my_ncombi_end,'my_nirred',my_nirred
2409 
2410   !COPY IRREDUCIBLE COMBINTATIONS TO BE DONE TO EACH PROCESSOR
2411   !my_nirred = my_ncombi_end - my_ncombi_start + 1
2412   ABI_MALLOC(my_list_combination_tmp,(power_disps(2),my_nirred))
2413   if(my_nirred /= 0)my_list_combination_tmp(:,:) = list_combination_tmp(:,my_ncombi_start:my_ncombi_end)
2414   ABI_MALLOC(my_index_irredcomb,(my_nirred))
2415   !DEALLOCATE list_combination_tmp
2416   ABI_FREE(list_combination_tmp)
2417 
2418   !COUNT SYMMETRIC COMBINATIONS TO IRREDUCIBLE COMBINATIONS ON EACH PROCESSOR
2419 
2420 
2421  if (my_nirred /= 0 ) then
2422      my_ncombi = my_nirred*(nsym**(power_disps(2)))
2423  end if
2424 
2425   !Copy irreducible combinations from list_combination tmp into my_list_combination on rank i
2426   !Make my_list_combination large enough for all symmetric combinations
2427 
2428   !ABI_MALLOC(my_list_combination,(power_disps(2),my_ncombi))
2429   !my_list_combination = 0
2430 
2431 
2432 
2433   !ABI_MALLOC(my_index_irredcomb,(my_ncombi_end-my_ncombi_start+1))
2434 !  if(my_ncombi /= 0)then
2435 !    do i=1,my_nirred
2436 !       my_list_combination(:,1+(i-1)*(nsym**(power_disps(2)))) = my_list_combination_tmp(:,i)
2437 !       my_index_irredcomb(i) = 1+(i-1)*(nsym**(power_disps(2)))
2438 !    end do
2439 !  endif
2440 
2441   !FOREGET ABOUT MY_LIST_COMBINATION_TMP
2442   !COMPUTE SYMMETRIC COMBINATIONS
2443   if(my_ncombi /= 0)then
2444     do i=1,my_nirred
2445       associate(comb => my_list_combination_tmp(:, i))
2446         !ABI_MALLOC(dummylist,(0))
2447         ! Get number of strain and displacements for this term
2448         ndisp = 0
2449         nstrain = 0
2450         do ii = 1,power_disps(2)
2451           if(comb(ii) > 0 .and. comb(ii) <= ncoeff_symsym)then
2452             ndisp = ndisp + 1
2453           else if(comb(ii) >= ncoeff_symsym)then
2454             nstrain = nstrain + 1
2455           endif
2456         enddo
2457        compute_sym = .true.
2458        iterm = my_index_irredcomb(i)-1
2459        call computeSymmetricCombinations(my_array_combination,list_symcoeff,list_symstr,ndisp,nsym,&
2460 &                                        comb(:ndisp+nstrain),power_disps(2),&
2461 &                                        my_ncombi,ncoeff_symsym,nstr_sym,nstrain, &
2462 &                                        compatibleCoeffs,compute_sym,comm,only_even=need_only_even_power)
2463 
2464        !ABI_FREE(dummylist)
2465      end associate
2466     enddo
2467   endif
2468 
2469   ABI_FREE(my_list_combination_tmp)
2470 
2471   ! Delete double combinations on each processor
2472   if(need_verbose .and. my_nirred /= 0)then
2473      write(message,'(1a,I4)')' Reduce reducible symmetric combinations on processor: ', my_rank+1
2474      call wrtout(std_out,message,'PERS')
2475   endif
2476   ! TODO: put my_array_combination into my_list_combination.
2477   !call reduce_zero_combinations(my_list_combination)
2478   call my_array_combination%tostatic(my_list_combination, size1=power_disps(2))
2479   call my_array_combination%finalize()
2480 
2481 
2482   ! Gather the Results into list_combination_tmp
2483   my_ncombi = size(my_list_combination,2)
2484   ABI_MALLOC(irank_ncombi,(nproc))
2485   call xmpi_allgather(my_ncombi,irank_ncombi,comm,ierr)
2486   if(need_verbose)then
2487     write(message,'(1a)')' Reduction on all processors finished. Gather results.'
2488     call wrtout(std_out,message,'COLL')
2489   endif
2490   ABI_MALLOC(list_combination_tmp,(power_disps(2),sum(irank_ncombi)))
2491   ABI_MALLOC(offsets,(nproc))
2492   offsets(1) = 0
2493   do i=1,nproc
2494     offsets(i) = sum(irank_ncombi(:i-1))*power_disps(2)
2495   enddo
2496 
2497   list_combination_tmp = 0
2498   ABI_MALLOC(buffsize,(nproc))
2499   do i = 1,nproc
2500     buffsize(i) = irank_ncombi(i)*power_disps(2)
2501   enddo
2502   call xmpi_gatherv(my_list_combination,size(my_list_combination),list_combination_tmp,buffsize,offsets,master,comm,ierr)
2503 
2504   !Deallocation of variables inside need_symmetric
2505   ABI_FREE(buffsize)
2506   ABI_FREE(my_list_combination)
2507   ABI_FREE(my_index_irredcomb)
2508   ABI_FREE(irank_ncombi)
2509   ABI_FREE(offsets)
2510 endif !compute_symmetric
2511 
2512 !Deallocation of arrays outside need_symmetric
2513 ABI_FREE(compatibleCoeffs)
2514 
2515 if(iam_master)then
2516   call reduce_zero_combinations(list_combination_tmp)
2517   ncombination = size(list_combination_tmp,2)
2518   !ABI_MALLOC(index_irred,(1))
2519   !index_irred = 1
2520   if(need_spcoupling)then !Check irreducibility of strain-phonon terms
2521    if(need_verbose)then
2522      write(message,'(1a)')' Reduce reducible Strain-Phonon combinations on master'
2523      call wrtout(std_out,message,'COLL')
2524    endif
2525 
2526    block
2527      type(IrreducibleCombinations_T) :: irred_combinations
2528      call irred_combinations%init()
2529      do i=1, ncombination
2530        if(any(list_combination_tmp(:,i) > ncoeff_symsym))then
2531          irreducible=irred_combinations%add_irr(list_combination_tmp(:,i), &
2532            & list_symcoeff, list_symstr, ncoeff_symsym, nsym, power_disps(2))
2533          if(.not. irreducible) list_combination_tmp(:,i) = 0
2534        endif
2535      end do
2536      call reduce_zero_combinations(list_combination_tmp)
2537      ncombination = size(list_combination_tmp,2)
2538      call irred_combinations%free()
2539    end block
2540  endif
2541 
2542 
2543  !    do i=2,ncombination
2544  !       if(any(list_combination_tmp(:,i) > ncoeff_symsym))then
2545  !         !irreducible = check_irreducibility(list_combination_tmp(:,i),list_combination_tmp(:,:i-1),list_symcoeff,&
2546  ! &       !                                   list_symstr,ncoeff_symsym,nsym,i-1,power_disps(2),index_irred)
2547  !         !if(.not. irreducible) list_combination_tmp(:,i) = 0
2548  !       endif
2549  !    enddo
2550  !    call reduce_zero_combinations(list_combination_tmp)
2551  !    ncombination = size(list_combination_tmp,2)
2552  ! ABI_FREE(index_irred)
2553 end if !iam_master
2554 if(need_verbose)then
2555   write(message,'(1x,I0,1a)') ncombination,' irreducible combinations generated '
2556   call wrtout(std_out,message,'COLL')
2557   write(message,'(1a)') ' Finished generating irreducible combinations'
2558   call wrtout(std_out,message,'COLL')
2559 endif
2560 
2561  !ABI_FREE(xcart)
2562  !ABI_FREE(xred)
2563 
2564 !MPI
2565  if(need_verbose .and. nproc > 1)then
2566    write(message,'(a,a)') ch10,' Redistribute the combinations over the CPU'
2567    call wrtout(std_out,message,'COLL')
2568  end if
2569 
2570  call xmpi_bcast(ncombination, master, comm, ierr)
2571 
2572  ncoeff_alone = mod(ncombination,nproc)
2573  my_ncoeff = int(aint(real(ncombination,sp)/(nproc)))
2574 
2575  if(my_rank >= (nproc-ncoeff_alone)) then
2576    my_ncoeff = my_ncoeff  + 1
2577  end if
2578 
2579 !Set the buffsize for mpi scatterv
2580  ABI_MALLOC(buffsize,(nproc))
2581  ABI_MALLOC(buffdispl,(nproc))
2582  do ii = 1,nproc
2583    buffsize(ii) = int(aint(real(ncombination,sp)/(nproc))*power_disps(2))
2584    if(ii > (nproc-ncoeff_alone)) then
2585      buffsize(ii) = buffsize(ii) + power_disps(2)
2586    end if
2587  end do
2588 
2589  buffdispl(1) = 0
2590  do ii = 2,nproc
2591    buffdispl(ii) = buffdispl(ii-1) + buffsize(ii-1)
2592  end do
2593 
2594  ABI_MALLOC(list_combination,(power_disps(2),my_ncoeff))
2595  list_combination = 0
2596 
2597  my_size = my_ncoeff*power_disps(2)
2598  call xmpi_scatterv(list_combination_tmp,buffsize,buffdispl,list_combination,my_size,master,&
2599 &                   comm,ierr)
2600 
2601  ABI_FREE(buffdispl)
2602  ABI_FREE(buffsize)
2603  ABI_FREE(list_combination_tmp)
2604 
2605 !write(std_out,*) "DEBUG shape list_symcoeff(:,:,:) before termsfromlist:", shape(list_symcoeff)
2606  !Compute the coefficients from the list of combination
2607  if(need_verbose .and. nproc > 1)then
2608    write(message,'(1a)')' Compute the coefficients'
2609    call wrtout(std_out,message,'COLL')
2610  end if
2611  ABI_MALLOC(coeffs_tmp,(my_ncoeff))
2612  nterm      = nsym
2613  ndisp_max  = power_disps(2)
2614  ncoeff_max = my_ncoeff
2615  do ii=1,my_ncoeff
2616    ABI_MALLOC(terms,(nterm))
2617    !write(std_out,*) "DEBUG list_combination(:,",ii,"): ", list_combination(:,ii)
2618    block
2619      logical :: reverse(ndisp_max)
2620      reverse=.False.
2621      call generateTermsFromList(cell,list_combination(:,ii),list_symcoeff,list_symstr,ncoeff_symsym,&
2622 &                             ndisp_max,nrpt,nstr_sym,nsym,nterm,terms, reverse=reverse)
2623      call polynomial_coeff_init(one,nterm,coeffs_tmp(ii),terms(1:nterm),check=.true.)
2624    end block
2625    !  Free the terms array
2626    do iterm=1,nterm
2627      call polynomial_term_free(terms(iterm))
2628    end do
2629    ABI_FREE(terms)
2630  end do
2631  ABI_FREE(cell)
2632 
2633 
2634  ABI_FREE(list_combination)
2635  ABI_FREE(list_symcoeff)
2636  ABI_FREE(list_symstr)
2637 
2638 !Final tranfert
2639 !1- Count the total number of coefficient
2640  ncoeff = 0
2641  do icoeff=1,ncoeff_max
2642    if (abs(coeffs_tmp(icoeff)%coefficient) >tol16) then
2643      ncoeff = ncoeff + 1
2644    end if
2645  end do
2646 
2647 !Get the total number of coefficients
2648 !ncoeff_max is the number of total coefficients before the symetries check
2649 !ncoeff_tot is the number of total coefficients after the symetries check
2650  ncoeff_tot = ncoeff!set the output
2651  call xmpi_sum(ncoeff_tot,comm,ierr)
2652  call xmpi_sum(ncoeff_max,comm,ierr)
2653 
2654 
2655 !Need to redistribute the coefficients over the CPU
2656 !Get the list with the number of coeff on each CPU
2657 !In order to be abble to compute the my_coeffindexes array which is for example:
2658 ! if CPU0 has 200  Coeff and CPU1 has 203 Coeff then
2659 ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403
2660  if(need_verbose .and. nproc > 1)then
2661    write(message,'(1a)')' Redistribute the coefficients over the CPU'
2662    call wrtout(std_out,message,'COLL')
2663  end if
2664 
2665  ABI_MALLOC(buffdispl,(nproc))
2666  buffdispl = 0
2667  buffdispl(my_rank+1) = my_ncoeff
2668  call xmpi_sum(buffdispl,comm,ierr)
2669  ABI_MALLOC(my_coeffindexes,(my_ncoeff))
2670  ABI_MALLOC(my_coefflist,(my_ncoeff))
2671  my_coeffindexes = 0
2672  my_coefflist = 0
2673  do icoeff=1,my_ncoeff
2674    my_coefflist(icoeff) = icoeff
2675    if(my_rank==0) then
2676      my_coeffindexes(icoeff) = icoeff
2677    else
2678      my_coeffindexes(icoeff) = sum(buffdispl(1:my_rank)) + icoeff
2679    end if
2680  end do
2681  ABI_FREE(buffdispl)
2682 
2683 !Compute the new number of coefficient per CPU
2684  if(need_distributed) then
2685    ncoeff_alone = mod(ncoeff_tot,nproc)
2686    my_newncoeff = int(aint(real(ncoeff_tot,sp)/(nproc)))
2687    if(my_rank >= (nproc-ncoeff_alone)) then
2688      my_newncoeff = my_newncoeff  + 1
2689    end if
2690  else
2691    my_newncoeff = ncoeff_tot
2692  end if
2693 
2694  ncoeff = my_newncoeff ! Set the output
2695 
2696 !2:compute the number of coefficients and the list of the corresponding
2697 !  coefficients for each CPU.
2698  ABI_MALLOC(my_newcoeffindexes,(my_newncoeff))
2699  if(need_distributed) then
2700    do icoeff=1,my_newncoeff
2701      if(my_rank >= (nproc-ncoeff_alone))then
2702        my_newcoeffindexes(icoeff)=int(aint(real(ncoeff_tot,sp)/(nproc)))*(my_rank)+&
2703 &                              (my_rank - (nproc-ncoeff_alone)) + icoeff
2704      else
2705        my_newcoeffindexes(icoeff)=(my_newncoeff)*(my_rank)  + icoeff
2706      end if
2707    end do
2708  else
2709    do icoeff=1,my_newncoeff
2710      my_newcoeffindexes(icoeff) = icoeff
2711    end do
2712  end if
2713 
2714 !2- Transfer
2715  if(.not.need_distributed)then
2716    if(.not.allocated(coefficients))then
2717      ABI_MALLOC(coefficients,(my_newncoeff))
2718    end if
2719  end if
2720  icoeff  = 0! icoeff is the current index in the total list of coefficients
2721  icoeff2 = 0! icoeff2 is the current index in the output coefficients array on each CPU
2722  icoeff3 = 0! icoeff3 is the current index in total new list of coefficients
2723  rank_to_send_save = 0
2724 
2725  do icoeff=1,ncoeff_max
2726 !  Need to send the rank with the chosen coefficient
2727    rank_to_send = 0
2728    my_icoeff = 0
2729    do ii=1,my_ncoeff
2730      if (my_coeffindexes(ii)==icoeff) then
2731        my_icoeff = ii
2732        if (abs(coeffs_tmp(my_icoeff)%coefficient) > tol16)then
2733          rank_to_send = my_rank
2734        else
2735          rank_to_send = -1
2736 !        Free the coefficient
2737          call polynomial_coeff_free(coeffs_tmp(ii))
2738        end if
2739        exit
2740      end if
2741    end do
2742    call xmpi_sum(rank_to_send, comm, ierr)
2743 !  This coefficient is not compute
2744    if (rank_to_send == -1) cycle
2745 
2746 !  increase icoeff3
2747    icoeff3 = icoeff3 + 1
2748 
2749 !  Find the receiver CPU
2750    rank_to_receive = 0
2751    do ii=1,my_newncoeff
2752      if (my_newcoeffindexes(ii)==icoeff3) then
2753        rank_to_receive = my_rank
2754      end if
2755    end do
2756    call xmpi_sum(rank_to_receive, comm, ierr)
2757 
2758    if(need_distributed.and.rank_to_send /= rank_to_send_save) then
2759      if(my_rank == rank_to_send_save)then
2760        ABI_FREE(coeffs_tmp)!Free memory if the current CPU has already distribute
2761                                           !all its own coefficients
2762      end if
2763      rank_to_send_save = rank_to_send
2764    end if
2765 
2766    if(need_distributed.and.my_rank == rank_to_receive)then
2767      if(.not.allocated(coefficients))then
2768        ABI_MALLOC(coefficients,(my_newncoeff))
2769      end if
2770    end if
2771 
2772 
2773    if (need_distributed)then
2774      if(my_rank==rank_to_send)then
2775        if(any(my_newcoeffindexes(:)==icoeff3))then
2776          icoeff2 = icoeff2 + 1
2777 !        Get the name of this coefficient
2778          call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2779          call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2780 &                                  coeffs_tmp(my_icoeff)%terms,name=name,check=.false.)
2781        else
2782          call polynomial_coeff_MPIsend(coeffs_tmp(my_icoeff), icoeff, rank_to_receive, comm)
2783        end if
2784 !      Free the coefficient
2785        call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2786      else
2787        if(any(my_newcoeffindexes(:)==icoeff3))then
2788          icoeff2 = icoeff2 + 1
2789          call polynomial_coeff_MPIrecv(coefficients(icoeff2), icoeff, rank_to_send, comm)
2790          call polynomial_coeff_getName(name,coefficients(icoeff2),symbols,recompute=.TRUE.)
2791          call polynomial_coeff_SetName(name,coefficients(icoeff2))
2792        end if
2793      end if
2794    else
2795      icoeff2 = icoeff2 + 1
2796 !    Get the name of this coefficient
2797      if(my_rank==rank_to_send)then
2798        call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2799        call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2800 &                                 coeffs_tmp(my_icoeff)%terms,name=name,check=.false.)
2801 !      Free the coefficient
2802        call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2803      end if
2804      call polynomial_coeff_broadcast(coefficients(icoeff2),rank_to_send, comm)
2805    end if
2806  end do
2807 
2808 !Debug write xml
2809 !     write (filename, "(A9,I2,A4)") "terms_set", my_rank+1,".xml"
2810 !     call polynomial_coeff_writeXML(coefficients,my_newncoeff,filename=filename)
2811 
2812  if(need_verbose)then
2813    write(message,'(1x,I0,2a)') ncoeff_tot,' coefficients generated ',ch10
2814    call wrtout(ab_out,message,'COLL')
2815    call wrtout(std_out,message,'COLL')
2816  end if
2817 
2818 
2819 !Final deallocation
2820  ABI_FREE(symbols)
2821  ABI_FREE(my_coeffindexes)
2822  ABI_FREE(my_newcoeffindexes)
2823  ABI_FREE(my_coefflist)
2824    ABI_SFREE(coeffs_tmp)
2825 
2826 end subroutine polynomial_coeff_getNorder

m_polynomial_coeff/polynomial_coeff_getOrder1 [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_getOrder1

FUNCTION

 Compute the first order polynomial coefficients from the list

INPUTS

 cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...)
 cutoff_in = cut-off for the inter atomic forces constants
 list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients,
                                    for each coefficients (ncoeff_sym), we store the symmetrics(nsym)
                                    the 6th first dimensions are :
                                       1 = direction of the IFC
                                       2 = index of the atom number 1 (1=>natom)
                                       3 = index of the atom number 2 (1=>natom)
                                       4 = indexes of the cell of the second atom
                                           (the atom number 1 is always in the cell 0 0 0)
                                       5 = weight of the term (-1 or 1)
                                       6 = indexes of the symmetric
 natom = number of atoms in the unit cell
 nrpt = number of cell
 nsym = number of symmetries in the system
 symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...)
 comm = MPI communicator

OUTPUT

 polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with
                                                              the polynomial_coeff
 ncoeff_out = number of coefficients

SOURCE

3754 subroutine polynomial_coeff_getOrder1(cell,coeffs_out,list_symcoeff,&
3755 &                                     natom,ncoeff_out,ncoeff,nrpt,nsym,&
3756 &                                     symbols)
3757 
3758  implicit none
3759 
3760 !Arguments ------------------------------------
3761 !scalars
3762  integer,intent(in)  :: natom,ncoeff,nsym,nrpt
3763  integer,intent(out) :: ncoeff_out
3764 !arrays
3765  integer,intent(in) :: cell(3,nrpt)
3766  integer,intent(in) :: list_symcoeff(6,ncoeff,nsym)
3767  character(len=5),intent(in) :: symbols(natom)
3768  type(polynomial_coeff_type),allocatable,intent(inout) :: coeffs_out(:)
3769 !Local variables-------------------------------
3770 !scalar
3771  integer :: ia,ib,icoeff,icoeff_tmp,irpt,irpt_ref
3772  integer :: isym,iterm,mu,ncoeff_max,ndisp,nstrain,nterm_max
3773  real(dp):: coefficient,weight
3774 !arrays
3775  integer,allocatable :: atindx(:,:),cells(:,:,:),dir_int(:)
3776  integer,allocatable :: power_disps(:),power_strain(:),strain(:)
3777  character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
3778  character(len=200):: name
3779  character(len=500) :: message
3780  type(polynomial_term_type),dimension(:),allocatable :: terms
3781  type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
3782 !TEST_AM
3783  character(len=fnlen) :: filename
3784 !TEST_AM
3785 ! *************************************************************************
3786 
3787 !Initialisation of variables
3788  nterm_max  = nsym
3789  ncoeff_max = ncoeff
3790  ndisp = 1
3791  nstrain = 0
3792  ABI_MALLOC(coeffs_tmp,(ncoeff_max))
3793  ABI_MALLOC(terms,(nterm_max))
3794 
3795 
3796  icoeff_tmp = 0
3797  ABI_MALLOC(atindx,(2,ndisp))
3798  ABI_MALLOC(cells,(3,2,ndisp))
3799  ABI_MALLOC(dir_int,(ndisp))
3800  ABI_MALLOC(power_disps,(ndisp))
3801  ABI_MALLOC(power_strain,(nstrain))
3802  ABI_MALLOC(strain,(nstrain))
3803 
3804 !Found the ref cell
3805  irpt_ref = 1
3806  do irpt=1,nrpt
3807    if(all(cell(:,irpt)==0))then
3808      irpt_ref = irpt
3809      exit
3810    end if
3811  end do
3812 
3813  write(message,'(3a)') " Irreductible coefficient and associated atom 1, atom 2 and direction:",ch10,&
3814 &                     " for the 1st order"
3815  call wrtout(std_out,message,'COLL')
3816 
3817  do icoeff=1,ncoeff
3818 !  Reset counter
3819    iterm = 0
3820    coefficient = one
3821    do isym=1,nsym
3822      ndisp   = 1
3823      nstrain = 0
3824      mu   = list_symcoeff(1,icoeff,isym)
3825      ia   = list_symcoeff(2,icoeff,isym)
3826      ib   = list_symcoeff(3,icoeff,isym)
3827      irpt = list_symcoeff(4,icoeff,isym)
3828      weight = list_symcoeff(5,icoeff,isym)
3829 !    Fill First term arrays
3830      atindx(1,1) = ia; atindx(2,1) = ib;
3831      dir_int(1) = mu
3832      power_disps(1)   = 1
3833      cells(:,1,1) = (/0,0,0/)
3834      cells(:,2,1) = cell(:,irpt)
3835      iterm = iterm + 1
3836      call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(iterm),&
3837 &                              power_disps,power_strain,strain,weight,check=.true.)
3838    end do!end do sym
3839 
3840    if(iterm > 0)then
3841 !  increase coefficients and set it
3842      icoeff_tmp = icoeff_tmp + 1
3843      call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),terms(1:iterm),check=.true.)
3844    end if
3845 
3846 !  Deallocate the terms
3847    do iterm=1,nterm_max
3848      call polynomial_term_free(terms(iterm))
3849    end do
3850  end do!end do coeff_sym
3851 
3852  ABI_FREE(terms)
3853  ABI_FREE(atindx)
3854  ABI_FREE(cells)
3855  ABI_FREE(dir_int)
3856  ABI_FREE(power_disps)
3857  ABI_FREE(power_strain)
3858  ABI_FREE(strain)
3859 
3860 !Count the number of terms
3861  ncoeff_out = 0
3862  do icoeff_tmp=1,ncoeff_max
3863    if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3864      ncoeff_out = ncoeff_out + 1
3865    end if
3866  end do
3867 
3868 !Transfer in the final array
3869  ABI_MALLOC(coeffs_out,(ncoeff_out))
3870  icoeff = 0
3871  do icoeff_tmp=1,ncoeff_max
3872    if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3873 !    Get the name of this coefficient
3874      call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
3875 !    Increase icoeff and fill the coeffs_out array
3876      icoeff = icoeff + 1
3877       call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
3878  &                               coeffs_out(icoeff),coeffs_tmp(icoeff_tmp)%terms,name=name)
3879 
3880      write(message,'(2a)')' ',trim(name)
3881      call wrtout(std_out,message,'COLL')
3882 
3883      do iterm = 1,coeffs_tmp(icoeff_tmp)%nterm
3884        write(message,'(a,I0,a,I0,2a)') '    Atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(1,1),&
3885 &                       ' and atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(2,1),&
3886 &                       ' in the direction ',mutodir(coeffs_tmp(icoeff_tmp)%terms(iterm)%direction(1))
3887        if(any(coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(:,2,1)/=0))then
3888          write(message,'(2a,I0,a,I0,a,I0,a)') trim(message),' in the cell ',&
3889 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(1,2,1),' ',&
3890 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(2,2,1),' ',&
3891 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(3,2,1),'.'
3892        end if
3893        call wrtout(std_out,message,'COLL')
3894      end do
3895    end if
3896  end do
3897 
3898 !TEST_AM
3899  filename = "terms_1st_order.xml"
3900  call polynomial_coeff_writeXML(coeffs_out,ncoeff_out,filename=filename)
3901 !TEST_AM
3902 
3903  write(message,'(a,1x,I0,a)') ch10,&
3904 &       ncoeff_out,' fitted coefficients for the 1st order '
3905  call wrtout(ab_out,message,'COLL')
3906  call wrtout(std_out,message,'COLL')
3907 
3908 !Deallocation
3909  do icoeff=1,ncoeff_max
3910    call polynomial_coeff_free(coeffs_tmp(icoeff))
3911  end do
3912  ABI_FREE(coeffs_tmp)
3913 
3914 end subroutine polynomial_coeff_getOrder1

m_polynomial_coeff/polynomial_coeff_init [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_init

FUNCTION

 Initialize a polynomial_coeff datatype

INPUTS

  name     = Name of the polynomial_coeff (Sr_y-O1_y)^3) for example
  nterm    = Number of terms (short range interaction) for this polynomial_coeff
  coefficient  = Value of the coefficient of this term
  terms(nterm)<type(polynomial_term)> = array of polynomial_term_type
  check   = TRUE if this list of terms has to be check. We remove the symetric of equivalent terms
                  for example:  ((Sr_y-O1_y)^1 and -1*(Sr_y-O1_y)^1 => zero
            FALSE, defaut, do nothing

OUTPUT

   polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype to be initialized

SOURCE

196 subroutine polynomial_coeff_init(coefficient,nterm,polynomial_coeff,terms,name,check)
197 
198  implicit none
199 
200 !Arguments ------------------------------------
201 !scalars
202  integer, intent(in) :: nterm
203  real(dp),intent(in) :: coefficient
204  logical,optional,intent(in) :: check
205 !arrays
206  character(len=200),optional,intent(in) :: name
207  type(polynomial_term_type),intent(in) :: terms(nterm)
208  type(polynomial_coeff_type), intent(out) :: polynomial_coeff
209 !Local variables-------------------------------
210 !scalar
211  integer :: iterm1,iterm2
212  integer :: ii,nterm_tmp
213  real(dp):: coefficient_tmp
214  logical :: check_in
215 !arrays
216  real(dp) :: weights(nterm)
217  character(len=200) :: name_tmp
218 ! *************************************************************************
219 !First free before initilisation
220  call polynomial_coeff_free(polynomial_coeff)
221  check_in = .false.
222  if(present(check)) check_in = check
223  if(check_in)then
224 !  Check if the list of term is available or contains identical terms
225 !  in this case, remove all the not needed terms
226    nterm_tmp = 0
227    weights(:) = one
228    do iterm1=1,nterm
229      if(abs(weights(iterm1)) < tol16)cycle ! FIXME: do nothing?
230      weights(iterm1) = terms(iterm1)%weight
231      do iterm2=iterm1+1,nterm
232        if(abs(weights(iterm2)) < tol16)cycle
233 !      if the terms are identical we check the weight
234        if(terms(iterm1)==terms(iterm2))then
235          weights(iterm1) = weights(iterm1) + terms(iterm2)%weight
236          weights(iterm2) = 0
237        end if
238      end do
239      if(abs(weights(iterm1)) > tol16) then
240        ! FIXME: weights(iterm1) = 1 ?
241        weights(iterm1)= anint(weights(iterm1)/weights(iterm1))
242      end if
243    end do
244 
245 
246 !  Count the number of terms
247    ! TODO: replace by: ??
248    !nterm_tmp=count(abs(weights) > tol16)
249    nterm_tmp = 0
250    do iterm1=1,nterm
251      if(abs(weights(iterm1)) > tol16)then
252        nterm_tmp = nterm_tmp + 1
253      end if
254    end do
255 
256    if (nterm_tmp ==0)then
257      coefficient_tmp = 0.0
258    else
259      coefficient_tmp = coefficient
260    end if
261 
262  else
263    nterm_tmp = nterm
264    coefficient_tmp = coefficient
265    weights(:) = terms(:)%weight
266  end if!end Check
267 
268  if(present(name))then
269    name_tmp = name
270  else
271    name_tmp = ""
272  end if
273 
274 
275 !Initilisation
276  polynomial_coeff%name = name_tmp
277  polynomial_coeff%nterm = nterm_tmp
278  polynomial_coeff%coefficient = coefficient_tmp
279  ABI_MALLOC(polynomial_coeff%terms,(polynomial_coeff%nterm))
280  iterm1 = 0
281  do ii = 1,nterm
282    if(abs(weights(ii)) > tol16)then
283      iterm1 = iterm1 + 1
284      call polynomial_term_init(terms(ii)%atindx,terms(ii)%cell,terms(ii)%direction,terms(ii)%ndisp,&
285 &                              terms(ii)%nstrain,polynomial_coeff%terms(iterm1),terms(ii)%power_disp,&
286 &                              terms(ii)%power_strain,terms(ii)%strain,terms(ii)%weight)
287    end if
288  end do
289 end subroutine polynomial_coeff_init

m_polynomial_coeff/polynomial_coeff_list_free [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_list_free

FUNCTION

 Free polynomial_coeff datatype

INPUTS

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

OUTPUT

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

SOURCE

353 subroutine polynomial_coeff_list_free(polynomial_coeff_list)
354 
355  implicit none
356 
357 !Arguments ------------------------------------
358 !scalars
359 !arrays
360  type(polynomial_coeff_type),allocatable, intent(inout) :: polynomial_coeff_list(:)
361 !Local variables-------------------------------
362 !scalar
363  integer :: i,ncoeff
364 !arrays
365 
366 ! *************************************************************************
367 
368 !Free output
369  if(allocated(polynomial_coeff_list))then
370    ncoeff = size(polynomial_coeff_list)
371    do i=1,ncoeff
372       call polynomial_coeff_free(polynomial_coeff_list(i))
373     enddo
374    ABI_FREE(polynomial_coeff_list)
375  endif
376 
377 end subroutine polynomial_coeff_list_free

m_polynomial_coeff/polynomial_coeff_MPIrecv [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_MPIrecv

FUNCTION

  MPI receive the polynomial_coefficent datatype

INPUTS

   tag = tag of the message to receive
   source = rank of Source
   comm = MPI communicator

SIDE EFFECTS

   coefficients<type(polynomial_coefficent_type)>=  polynomial_coeff datatype

SOURCE

755 subroutine polynomial_coeff_MPIrecv(coefficients, tag, source, comm)
756 
757  implicit none
758 
759 !Arguments ------------------------------------
760 !array
761  type(polynomial_coeff_type),intent(inout) :: coefficients
762  integer, intent(in) :: source,comm,tag
763 
764 !Local variables-------------------------------
765 !scalars
766  integer :: ierr,ii
767 !arrays
768 
769 ! *************************************************************************
770 
771  if (xmpi_comm_size(comm) == 1) return
772 
773 
774 ! Free the output
775   call polynomial_coeff_free(coefficients)
776 
777  ! Transmit variables
778   call xmpi_recv(coefficients%name, source, 9*tag+0, comm, ierr)
779   call xmpi_recv(coefficients%nterm, source, 9*tag+1, comm, ierr)
780   call xmpi_recv(coefficients%coefficient, source, 9*tag+2, comm, ierr)
781 
782  !Allocate arrays on the other nodes.
783   ABI_MALLOC(coefficients%terms,(coefficients%nterm))
784   do ii=1,coefficients%nterm
785     call polynomial_term_free(coefficients%terms(ii))
786   end do
787 
788 ! Set the number of term on each node (needed for allocations of array)
789   do ii = 1,coefficients%nterm
790     call xmpi_recv(coefficients%terms(ii)%ndisp, source, 9*tag+3, comm, ierr)
791     call xmpi_recv(coefficients%terms(ii)%nstrain, source, 9*tag+4, comm, ierr)
792   end do
793 
794 ! Allocate arrays on the other nodes
795   do ii = 1,coefficients%nterm
796     ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
797     coefficients%terms(ii)%atindx = 0
798     ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
799     ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
800     ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
801     ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
802     ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
803   end do
804 
805 ! Transfert value
806   do ii = 1,coefficients%nterm
807     call xmpi_recv(coefficients%terms(ii)%weight, source, 9*tag+5, comm, ierr)
808     call xmpi_recv(coefficients%terms(ii)%atindx, source, 9*tag+6, comm, ierr)
809     call xmpi_recv(coefficients%terms(ii)%direction, source, 9*tag+7, comm, ierr)
810     call xmpi_recv(coefficients%terms(ii)%cell, source, 9*tag+8, comm, ierr)
811     call xmpi_recv(coefficients%terms(ii)%power_disp, source, 9*tag+9, comm, ierr)
812     call xmpi_recv(coefficients%terms(ii)%power_strain, source, 9*tag+10, comm, ierr)
813     call xmpi_recv(coefficients%terms(ii)%strain, source, 9*tag+11, comm, ierr)
814   end do
815 
816 end subroutine polynomial_coeff_MPIrecv

m_polynomial_coeff/polynomial_coeff_MPIsend [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_MPIsend

FUNCTION

  MPI send the polynomial_coefficent datatype

INPUTS

   tag = tag of the message to send
   dest= rank of Dest
   comm= MPI communicator

SIDE EFFECTS

   polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

SOURCE

693 subroutine polynomial_coeff_MPIsend(coefficients, tag, dest, comm)
694 
695  implicit none
696 
697 !Arguments ------------------------------------
698 !array
699  type(polynomial_coeff_type),intent(inout) :: coefficients
700  integer, intent(in) :: dest,comm,tag
701 
702 !Local variables-------------------------------
703 !scalars
704  integer :: ierr,ii
705  integer :: my_rank
706 !arrays
707 
708 ! *************************************************************************
709 
710  if (xmpi_comm_size(comm) == 1) return
711 
712   my_rank = xmpi_comm_rank(comm)
713 ! Transmit variables
714   call xmpi_send(coefficients%name, dest, 9*tag+0, comm, ierr)
715   call xmpi_send(coefficients%nterm, dest, 9*tag+1, comm, ierr)
716   call xmpi_send(coefficients%coefficient, dest, 9*tag+2, comm, ierr)
717 
718 ! Set the number of term on each node (needed for allocations of array)
719   do ii = 1,coefficients%nterm
720     call xmpi_send(coefficients%terms(ii)%ndisp, dest, 9*tag+3, comm, ierr)
721     call xmpi_send(coefficients%terms(ii)%nstrain, dest, 9*tag+4, comm, ierr)
722   end do
723 
724 ! Transfert value
725   do ii = 1,coefficients%nterm
726       call xmpi_send(coefficients%terms(ii)%weight, dest, 9*tag+5, comm, ierr)
727       call xmpi_send(coefficients%terms(ii)%atindx, dest, 9*tag+6, comm, ierr)
728       call xmpi_send(coefficients%terms(ii)%direction, dest, 9*tag+7, comm, ierr)
729       call xmpi_send(coefficients%terms(ii)%cell, dest, 9*tag+8, comm, ierr)
730       call xmpi_send(coefficients%terms(ii)%power_disp, dest, 9*tag+9, comm, ierr)
731       call xmpi_send(coefficients%terms(ii)%power_strain, dest, 9*tag+10, comm, ierr)
732       call xmpi_send(coefficients%terms(ii)%strain, dest, 9*tag+11, comm, ierr)
733   end do
734 
735 end subroutine polynomial_coeff_MPIsend

m_polynomial_coeff/polynomial_coeff_setCoefficient [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_setCoefficient

FUNCTION

 set the coefficient for of polynomial_coeff

INPUTS

 coefficient = coefficient of this coefficient

OUTPUT

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

SOURCE

398 subroutine polynomial_coeff_setCoefficient(coefficient,polynomial_coeff)
399 
400  implicit none
401 
402 !Arguments ------------------------------------
403 !scalars
404  real(dp),intent(in) :: coefficient
405 !arrays
406  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
407 !Local variables-------------------------------
408 !scalar
409 !arrays
410 ! *************************************************************************
411 
412  polynomial_coeff%coefficient = coefficient
413 
414 end subroutine polynomial_coeff_setCoefficient

m_polynomial_coeff/polynomial_coeff_setName [ Functions ]

[ Top ] [ Functions ]

NAME

 polynomial_coeff_setName

FUNCTION

 set the name of a  polynomial_coeff type

INPUTS

 name = name of the coeff

OUTPUT

 polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype

SOURCE

433 subroutine polynomial_coeff_setName(name,polynomial_coeff)
434 
435  implicit none
436 
437 !Arguments ------------------------------------
438 !scalars
439 !arrays
440  character(len=200),intent(in) :: name
441  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
442 !Local variables-------------------------------
443 !scalar
444 !arrays
445 ! *************************************************************************
446 
447  polynomial_coeff%name = name
448 
449 end subroutine polynomial_coeff_setName

m_polynomial_coeff/polynomial_coeff_type [ Types ]

[ Top ] [ Types ]

NAME

 polynomial_coeff_type

FUNCTION

 structure for a polynomial coefficient
 contains the value of the coefficient and a
 list of terms (displacements and/or strain) relating to the coefficient by symmetry

SOURCE

133  type, public :: polynomial_coeff_type
134 
135    character(len=200) :: name = ""
136 !     Name of the polynomial_coeff (Sr_y-O1_y)^3) for example
137 
138    integer :: nterm = 0
139 !     Number of terms (short range interaction) for this polynomial_coeff
140 
141    real(dp) :: coefficient = zero
142 !     coefficient = value of the coefficient of this term
143 !     \frac{\partial E^{k}}{\partial \tau^{k}}
144 
145    type(polynomial_term_type),dimension(:),allocatable :: terms
146 !     polynomial_term(nterm)<type(polynomial_term)>
147 !     contains all the displacements for this coefficient
148 
149  end type polynomial_coeff_type

m_polynomial_coeff/reduce_zero_combinations [ Functions ]

[ Top ] [ Functions ]

NAME

 reduce_zero_combinations

FUNCTION

 Sort out list of zeros in a list of integers

INPUTS

OUTPUT

SOURCE

4472 subroutine reduce_zero_combinations(combination_list)
4473 !Arguments ------------------------------------
4474  implicit none
4475 
4476 !Arguments ------------------------------------
4477  !scalar
4478  !array
4479  integer,allocatable,intent(inout) :: combination_list(:,:)
4480 !local
4481 !variable
4482   integer :: i,j
4483   integer,allocatable :: combination_list_tmp(:,:)
4484 !array
4485 ! *************************************************************************
4486 
4487  !Reduce zero combinations
4488  i = 0
4489  do j=1,size(combination_list,2)
4490    if(any(combination_list(:,j) /= 0))then
4491       i = i + 1
4492    endif
4493  enddo
4494 
4495  ABI_MALLOC(combination_list_tmp,(size(combination_list,1),i))
4496  i = 0
4497  do j=1,size(combination_list,2)
4498    if(any(combination_list(:,j) /= 0))then
4499       i = i + 1
4500       combination_list_tmp(:,i) = combination_list(:,j)
4501    endif
4502  enddo
4503  ABI_FREE(combination_list)
4504  ABI_MALLOC(combination_list,(size(combination_list_tmp,1),i))
4505  combination_list = combination_list_tmp
4506  ABI_FREE(combination_list_tmp)
4507 
4508 end subroutine reduce_zero_combinations

m_polynomial_coeff/sort_combination [ Functions ]

[ Top ] [ Functions ]

NAME

 sort_combination

FUNCTION

 Sort a list of integer from small to large if it contains zeros will be put to highest indexe

INPUTS

OUTPUT

SOURCE

4260 subroutine sort_combination(combination,n_int)
4261 !Arguments ------------------------------------
4262  implicit none
4263 
4264 !Arguments ------------------------------------
4265  !scalar
4266  integer,intent(in) :: n_int
4267  !array
4268  integer,intent(inout) :: combination(n_int)
4269 !local
4270 !variable
4271   integer :: j,k,cnt,tmp_int1,tmp_int2
4272 !array
4273 ! *************************************************************************
4274 
4275 j=2
4276 do while(j <= n_int)
4277    k = j
4278    cnt = 1
4279    do while(k >= 2 .and. cnt == 1)
4280       if(combination(k-1) > combination(k) .and. combination(k) > 0)then
4281            tmp_int1 = combination(k-1)
4282            tmp_int2 = combination(k)
4283            combination(k) = tmp_int1
4284            combination(k-1) = tmp_int2
4285            k=k-1
4286       else
4287         cnt = cnt + 1
4288       end if
4289    end do
4290    j = j+1
4291 end do
4292 
4293 end subroutine sort_combination

m_polynomial_coeff/sort_combination_list [ Functions ]

[ Top ] [ Functions ]

NAME

 sort_combination_list

FUNCTION

 Sort a list of integer list from small to large if it contains zeros will be put to highest index

INPUTS

OUTPUT

SOURCE

4310 subroutine sort_combination_list(combination_list,n_int,n_list)
4311 !Arguments ------------------------------------
4312  implicit none
4313 
4314 !Arguments ------------------------------------
4315  !scalar
4316  integer,intent(in) :: n_int,n_list
4317  !array
4318  integer,intent(inout) :: combination_list(n_int,n_list)
4319 !local
4320 !variable
4321   integer :: i
4322 !array
4323 ! *************************************************************************
4324 
4325 do i=1,n_list
4326    call sort_combination(combination_list(:,i),n_int)
4327 end do
4328 
4329 end subroutine sort_combination_list