TABLE OF CONTENTS


ABINIT/m_polynomial_coeff [ Functions ]

[ Top ] [ Functions ]

NAME

 m_polynomial_coeff

FUNCTION

 Module with the datatype polynomial coefficients

COPYRIGHT

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

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_polynomial_coeff
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_polynomial_term
32  use m_xmpi
33 #ifdef HAVE_MPI2
34  use mpi
35 #endif
36 
37  use m_sort,      only : sort_dp
38  use m_io_tools,  only : open_file, get_unit
39  use m_symtk,     only : symchk, symatm
40  use m_crystal,   only : crystal_t,symbols_crystal
41  use m_supercell, only : getPBCIndexes_supercell,distance_supercell,findBound_supercell
42  use m_geometry,  only : xcart2xred,metric
43  use m_dtfil,     only : isfile
44 
45  implicit none
46 
47  public :: polynomial_coeff_broadcast
48  public :: polynomial_coeff_evaluate
49  public :: polynomial_coeff_free
50  public :: polynomial_coeff_init
51  public :: polynomial_coeff_getList
52  public :: polynomial_coeff_getName
53  public :: polynomial_coeff_getNorder
54  public :: polynomial_coeff_getOrder1
55  public :: polynomial_coeff_MPIrecv
56  public :: polynomial_coeff_MPIsend
57  public :: polynomial_coeff_setName
58  public :: polynomial_coeff_setCoefficient
59  public :: polynomial_coeff_writeXML
60  private :: computeNorder
61  private :: computeCombinationFromList
62  private :: getCoeffFromList
63  private :: generateTermsFromList

m_polynomial_coeff/coeffs_compare [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ Functions ]

NAME

  equal

FUNCTION

INPUTS

OUTPUT

SOURCE

3306 pure function coeffs_compare(c1,c2) result (res)
3307 !Arguments ------------------------------------
3308 
3309 !This section has been created automatically by the script Abilint (TD).
3310 !Do not modify the following lines by hand.
3311 #undef ABI_FUNC
3312 #define ABI_FUNC 'coeffs_compare'
3313 !End of the abilint section
3314 
3315  implicit none
3316 
3317 !Arguments ------------------------------------
3318   type(polynomial_coeff_type), intent(in) :: c1,c2
3319   logical :: res
3320 !local
3321 !variable
3322   integer :: iterm1,iterm2
3323 !array
3324   integer :: blkval(2,c1%nterm)
3325 ! *************************************************************************
3326   res = .false.
3327   blkval = 0
3328   do iterm1=1,c1%nterm
3329     if(blkval(1,iterm1)==1)cycle!already found
3330     do iterm2=1,c2%nterm
3331       if(blkval(2,iterm2)==1)cycle!already found
3332       if(c1%terms(iterm1)==c2%terms(iterm2)) then
3333         blkval(1,iterm1) = 1
3334         blkval(2,iterm2) = 1
3335       end if
3336     end do
3337   end do
3338   if(.not.any(blkval(:,:)==0))res = .true.
3339 
3340 end function coeffs_compare

m_polynomial_coeff/computeCombinationFromList [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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 we 1)
 max_power_strain = maximum order of the strain of the strain phonon coupling
 nmodel_tot = current number of combination already computed (start we 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 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
 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

PARENTS

CHILDREN

SOURCE

2741 recursive subroutine computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
2742 &                                  index_coeff_in,list_combination,icoeff,max_power_strain,nmodel_tot,&
2743 &                                  natom,ncoeff,nstr,nmodel,nrpt,nsym,power_disp,power_disp_min,&
2744 &                                  power_disp_max,symbols,nbody,only_odd_power,only_even_power,&
2745 &                                  compute,anharmstr,spcoupling)
2746 
2747 
2748 !This section has been created automatically by the script Abilint (TD).
2749 !Do not modify the following lines by hand.
2750 #undef ABI_FUNC
2751 #define ABI_FUNC 'computeCombinationFromList'
2752 !End of the abilint section
2753 
2754  implicit none
2755 
2756 !Arguments ---------------------------------------------
2757 !scalar
2758  integer,intent(in) :: natom,ncoeff,power_disp,power_disp_min,power_disp_max
2759  integer,intent(in) :: max_power_strain,nmodel,nsym,nrpt,nstr
2760  integer,intent(inout) :: icoeff,nmodel_tot
2761  logical,optional,intent(in) :: compute,anharmstr,spcoupling
2762  integer,optional,intent(in) :: nbody
2763  logical,optional,intent(in) :: only_odd_power,only_even_power
2764 !arrays
2765  integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff+nstr,ncoeff+nstr)
2766  integer,intent(in) :: list_coeff(6,ncoeff,nsym),list_str(nstr,nsym,2)
2767  integer,intent(in) :: index_coeff_in(power_disp-1)
2768  integer,intent(out) :: list_combination(power_disp_max,nmodel)
2769  character(len=5),intent(in) :: symbols(natom)
2770 !Local variables ---------------------------------------
2771 !scalar
2772  integer :: icoeff1,icoeff2,nbody_in,ii,jj
2773  logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling
2774  logical :: need_only_odd_power,need_only_even_power
2775 !arrays
2776  integer :: powers(power_disp)
2777  integer,allocatable :: index_coeff(:)
2778 ! *************************************************************************
2779 
2780 !Set the inputs
2781  need_compute = .TRUE.
2782  need_anharmstr = .TRUE.
2783  need_spcoupling = .TRUE.
2784  need_only_odd_power = .FALSE.
2785  need_only_even_power = .FALSE.
2786  nbody_in = 0 !all kind of terms
2787  if(present(compute)) need_compute = compute
2788  if(present(nbody)) nbody_in = nbody
2789  if(present(anharmstr)) need_anharmstr = anharmstr
2790  if(present(spcoupling)) need_spcoupling = spcoupling
2791  if(present(only_odd_power)) need_only_odd_power = only_odd_power
2792  if(present(only_even_power)) need_only_even_power = only_even_power
2793 
2794  if(power_disp <= power_disp_max)then
2795 
2796 !  Initialisation of variables
2797    ABI_ALLOCATE(index_coeff,(power_disp))
2798    index_coeff(1:power_disp-1) = index_coeff_in(:)
2799 !  Loop over ncoeff+nstr
2800    do icoeff1=icoeff,ncoeff+nstr
2801 
2802 !    Reset the flag compatible and possible
2803      compatible = .TRUE.
2804      possible   = .TRUE.
2805 
2806 !    If the power_disp is one, we need to set icoeff to icoeff1
2807      if(power_disp==1) then
2808        icoeff = icoeff1
2809        if(compatibleCoeffs(icoeff,icoeff1)==0)then
2810          compatible = .FALSE.
2811        end if
2812      end if
2813 !    If the distance between the 2 coefficients is superior than the cut-off, we cycle.
2814     do icoeff2=1,power_disp-1
2815       if(compatibleCoeffs(index_coeff(icoeff2),icoeff1)==0)then
2816         compatible = .FALSE.
2817       end if
2818     end do
2819 
2820      if (.not.compatible) cycle !The distance is not compatible
2821 
2822 !    Set the index of the new coeff in the list
2823      index_coeff(power_disp) = icoeff1
2824 !    Do some checks
2825 !    -------------
2826 !    1-Check if the coefficient is full anharmonic strain and if we need to compute it
2827      if(all(index_coeff > ncoeff))then
2828        compatible = (need_anharmstr .or. need_spcoupling)
2829        possible = need_anharmstr
2830      end if
2831 !    2-Check if the coefficient is strain-coupling and if we need to compute it
2832      if(any(index_coeff < ncoeff) .and. any(index_coeff > ncoeff))then
2833        possible   = need_spcoupling
2834        compatible = need_spcoupling
2835        if(count(index_coeff > ncoeff) > max_power_strain)then
2836          possible = .false.
2837          compatible = .false.
2838        end if
2839      end if
2840 
2841 
2842      if(power_disp >= power_disp_min) then
2843 
2844 !      count the number of body
2845        powers(:) = 1
2846        do ii=1,power_disp
2847          do jj=ii+1,power_disp
2848            if (powers(jj) == 0) cycle
2849            if(index_coeff(ii)==index_coeff(jj))then
2850              powers(ii) = powers(ii) + 1
2851              powers(jj) = 0
2852            end if
2853          end do
2854        end do
2855 
2856 !      check the only_odd and only_even flags
2857        if(any(mod(powers(1:power_disp),2) /=0) .and. need_only_even_power) then
2858          possible = .false.
2859        end if
2860        if(any(mod(powers(1:power_disp),2) ==0) .and. need_only_odd_power)then
2861          possible = .false.
2862        end if
2863 !      Check the nbody flag
2864        if(nbody_in /= 0)then
2865          if(power_disp-count(powers==0) > nbody_in) then
2866            possible = .false.
2867            compatible = .false.
2868          end if
2869        end if
2870 
2871        if(possible) then
2872 !        increase coefficients and set it
2873          nmodel_tot = nmodel_tot + 1
2874          if(need_compute)then
2875            list_combination(1:power_disp,nmodel_tot) = index_coeff
2876          end if
2877        end if
2878      end if!end if power_disp < power_disp_min
2879 
2880 !    If the model is still compatbile with the input flags, we continue.
2881      if(compatible)then
2882        call computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,&
2883 &                                     index_coeff,list_combination,icoeff1,max_power_strain,&
2884 &                                     nmodel_tot,natom,ncoeff,nstr,nmodel,nrpt,nsym,power_disp+1,&
2885 &                                     power_disp_min,power_disp_max,symbols,nbody=nbody_in,&
2886 &                                     compute=need_compute,anharmstr=need_anharmstr,&
2887 &                                     spcoupling=need_spcoupling,only_odd_power=need_only_odd_power,&
2888 &                                     only_even_power=need_only_even_power)
2889      end if
2890    end do
2891    ABI_DEALLOCATE(index_coeff)
2892  end if
2893 
2894 end subroutine computeCombinationFromList

m_polynomial_coeff/computeNorder [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

CHILDREN

SOURCE

2489 recursive subroutine computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,&
2490 &                                  index_coeff_in,icoeff,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,&
2491 &                                  nrpt,nsym,power_disp,power_disp_min,power_disp_max,symbols,nbody,&
2492 &                                  compute,anharmstr,spcoupling,distributed)
2493 
2494 
2495 !This section has been created automatically by the script Abilint (TD).
2496 !Do not modify the following lines by hand.
2497 #undef ABI_FUNC
2498 #define ABI_FUNC 'computeNorder'
2499 !End of the abilint section
2500 
2501  implicit none
2502 
2503 !Arguments ---------------------------------------------
2504 !scalar
2505  integer,intent(in) :: natom,ncoeff,power_disp,power_disp_min,power_disp_max,ncoeff_out,nsym,nrpt,nstr
2506  integer,intent(inout) :: icoeff,icoeff_tot
2507  logical,optional,intent(in) :: compute,anharmstr,spcoupling,distributed
2508  integer,optional,intent(in) :: nbody
2509 !arrays
2510  integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff+nstr,ncoeff+nstr)
2511  integer,intent(in) :: list_coeff(6,ncoeff,nsym),list_str(nstr,nsym,2)
2512  integer,intent(in) :: index_coeff_in(power_disp-1)
2513  type(polynomial_coeff_type),intent(inout) :: coeffs_out(ncoeff_out)
2514  character(len=5),intent(in) :: symbols(natom)
2515 !Local variables ---------------------------------------
2516 !scalar
2517  integer :: ia,ib,ii,icoeff1,icoeff_tmp
2518  integer :: iterm,nbody_in,ncoeff_max,pa,pb
2519  integer :: ndisp_max,nterm_max
2520  real(dp):: coefficient
2521  logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling,need_distributed
2522 !arrays
2523  integer,allocatable :: index_coeff(:)
2524  character(len=200):: name
2525  type(polynomial_term_type),dimension(:),allocatable :: terms
2526  type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
2527 ! *************************************************************************
2528 
2529 !Set the inputs
2530  need_compute = .TRUE.
2531  need_anharmstr = .TRUE.
2532  need_spcoupling = .TRUE.
2533  need_distributed = .FALSE.
2534  nbody_in = 0 !all kind of terms
2535  if(present(compute)) need_compute = compute
2536  if(present(nbody)) nbody_in = nbody
2537  if(present(anharmstr)) need_anharmstr = anharmstr
2538  if(present(spcoupling)) need_spcoupling = spcoupling
2539  if(present(distributed)) need_distributed  = distributed
2540  if(power_disp <= power_disp_max)then
2541 
2542 !  Initialisation of variables
2543    nterm_max  = nsym
2544    ncoeff_max = (ncoeff+nstr)
2545    ndisp_max = power_disp
2546    icoeff_tmp = 0
2547    ABI_ALLOCATE(coeffs_tmp,(ncoeff_max))
2548    ABI_ALLOCATE(terms,(nterm_max))
2549    ABI_ALLOCATE(index_coeff,(power_disp))
2550 
2551    index_coeff(1:power_disp-1) = index_coeff_in(:)
2552 
2553    do icoeff1=icoeff,ncoeff+nstr
2554 !    If the distance between the 2 coefficients is superior than the cut-off,
2555 !    we cycle
2556 !    If the power_disp is one, we need to set icoeff to icoeff1
2557      if(power_disp==1) icoeff = icoeff1
2558 
2559      if(compatibleCoeffs(icoeff,icoeff1)==0) cycle
2560 
2561 !    Reset the flag compatible and possible
2562      compatible = .TRUE.
2563      possible   = .TRUE.
2564 
2565      index_coeff(power_disp) = icoeff1
2566      iterm = 0
2567      coefficient = one
2568 
2569      if(power_disp >= power_disp_min) then
2570        call generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,&
2571 &                                 ndisp_max,nrpt,nstr,nsym,iterm,terms)
2572 
2573        if(iterm > 0)then
2574 !        Do some checks
2575 !        -------------
2576 !        1-Check if the coefficient is full anharmonic strain and if we need to compute it
2577          if(terms(1)%ndisp == 0)then
2578            compatible = (need_anharmstr .or. need_spcoupling)
2579            possible = need_anharmstr
2580          end if
2581 !        1-Check if the coefficient is strain-coupling and if we need to compute it
2582          if(terms(1)%nstrain > 0.and.terms(1)%ndisp > 0)then
2583            possible   = need_spcoupling
2584            compatible = need_spcoupling
2585          end if
2586 
2587 !        ------------
2588 !        2-Check if this terms is compatible with nbody
2589          if(nbody_in > 0)then
2590            pa = 1 ; pb = 1
2591            ia = 0 ; ib = 0
2592 !          Count the number of terms and the power_disp
2593            do ii=1,terms(1)%ndisp
2594              if(terms(1)%nstrain > 0) then
2595                pb = pb*terms(1)%power_disp(ii)
2596                ib = ib + 1
2597              else
2598                pa = pa*terms(1)%power_disp(ii)
2599                ia = ia + 1
2600              end if
2601            end do
2602            if(ia <= nbody_in)then
2603              if(ia==nbody_in.and.abs(mod(pa,2)) < tol16)then
2604                if(ib==0)then
2605                  compatible = .FALSE.
2606                  possible   = .TRUE.
2607                else if (ib==nbody_in.and.abs(mod(pb,2)) < tol16) then
2608                  compatible = .FALSE.
2609                  possible   = .TRUE.
2610                else
2611                 possible = .FALSE.
2612                 compatible = .FALSE.
2613                end if
2614              else
2615                 possible = .FALSE.
2616                 compatible = .FALSE.
2617              end if
2618            else
2619              compatible = .FALSE.
2620              possible = .FALSE.
2621            end if
2622          end if
2623 
2624          if(possible)then
2625 !          increase coefficients and set it
2626            icoeff_tmp = icoeff_tmp + 1
2627            icoeff_tot = icoeff_tot + 1
2628            call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),&
2629 &                                     terms(1:iterm),check=.true.)
2630 
2631          end if
2632        end if
2633 
2634 !      Deallocate the terms
2635        do iterm=1,nterm_max
2636          call polynomial_term_free(terms(iterm))
2637        end do
2638 
2639      end if!end if power_disp < power_disp_min
2640 
2641      if(compatible)then
2642        call computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,index_coeff,&
2643 &                         icoeff1,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,nrpt,nsym,power_disp+1,&
2644 &                         power_disp_min,power_disp_max,symbols,nbody=nbody_in,compute=need_compute,&
2645 &                         anharmstr=need_anharmstr,spcoupling=need_spcoupling)
2646      end if
2647    end do
2648 
2649    ABI_DEALLOCATE(terms)
2650    ABI_DEALLOCATE(index_coeff)
2651 
2652 !  Transfer in the final array
2653    icoeff1 = 0
2654    do icoeff_tmp=1,ncoeff_max
2655      if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
2656 !      Increase icoeff and fill the coeffs_out array
2657        icoeff_tot = icoeff_tot + 1
2658        if(need_compute)then
2659          name = ''
2660 !        Get the name of this coefficient
2661          call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
2662          call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
2663 &                                   coeffs_out(icoeff_tot),coeffs_tmp(icoeff_tmp)%terms,&
2664 &                                   name=name)
2665        end if
2666      end if
2667    end do
2668 !  Deallocation
2669    do icoeff1=1,ncoeff_max
2670      call polynomial_coeff_free(coeffs_tmp(icoeff1))
2671    end do
2672    ABI_DEALLOCATE(coeffs_tmp)
2673  end if
2674 
2675 end subroutine computeNorder

m_polynomial_coeff/generateTermsFromList [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

3011 subroutine generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,ndisp_max,&
3012 &                                nrpt,nstr,nsym,nterm,terms)
3013 
3014 
3015 !This section has been created automatically by the script Abilint (TD).
3016 !Do not modify the following lines by hand.
3017 #undef ABI_FUNC
3018 #define ABI_FUNC 'generateTermsFromList'
3019 !End of the abilint section
3020 
3021  implicit none
3022 
3023 !Arguments ------------------------------------
3024 !scalar
3025  integer,intent(in) :: ndisp_max,ncoeff,nrpt,nstr,nsym
3026  integer,intent(out):: nterm
3027 !arrays
3028  integer,intent(in) :: index_coeff(ndisp_max)
3029  integer,intent(in) :: cell(3,nrpt),list_coeff(6,ncoeff,nsym)
3030  integer,intent(in) :: list_str(nstr,nsym,2)
3031  type(polynomial_term_type),intent(out) :: terms(nsym)
3032 !Local variables-------------------------------
3033 !scalar
3034  integer :: ia,ib,icoeff_str,idisp,irpt
3035  integer :: isym,ndisp,nstrain,mu
3036  real(dp):: weight
3037 !arrays
3038  integer :: atindx(2,ndisp_max),cells(3,2,ndisp_max),dir_int(ndisp_max),strain(ndisp_max)
3039  integer :: power_disps(ndisp_max),power_strain(ndisp_max)
3040 
3041 ! *************************************************************************
3042  nterm = 0
3043 !Loop over symetries
3044  do isym=1,nsym
3045 !Treat this coeff
3046    weight = 1
3047    ndisp = 0
3048    nstrain = 0
3049    do idisp=1,ndisp_max
3050 !    Get index of this displacement term
3051 !    Check if the index is not zero
3052      if(index_coeff(idisp)==0)cycle
3053      if(index_coeff(idisp)<=ncoeff)then
3054        ndisp = ndisp + 1
3055        mu   = list_coeff(1,index_coeff(idisp),isym)
3056        ia   = list_coeff(2,index_coeff(idisp),isym)
3057        ib   = list_coeff(3,index_coeff(idisp),isym)
3058        irpt = list_coeff(4,index_coeff(idisp),isym)
3059        weight = weight*list_coeff(5,index_coeff(idisp),isym)
3060 !      Fill First term arrays
3061        atindx(1,idisp) = ia; atindx(2,idisp) = ib;
3062        dir_int(idisp) = mu
3063        power_disps(idisp)   = 1
3064        cells(:,1,idisp) = (/0,0,0/)
3065        cells(:,2,idisp) = cell(:,irpt)
3066      else
3067        nstrain = nstrain + 1
3068        icoeff_str = index_coeff(idisp)-ncoeff
3069        strain(nstrain) = list_str(icoeff_str,isym,1)
3070        power_strain(nstrain)  = 1
3071        weight = weight*list_str(icoeff_str,isym,2)
3072      end if
3073    end do
3074    nterm = nterm + 1
3075    call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(nterm),power_disps,&
3076 &                            power_strain,strain,weight,check=.true.)
3077  end do!end do sym
3078 
3079 
3080 end subroutine generateTermsFromList

m_polynomial_coeff/getCoeffFromList [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_polynomial_coeff

CHILDREN

SOURCE

2932 function getCoeffFromList(list_coeff,ia,ib,irpt,mu,ncoeff) result(coeff)
2933 
2934 
2935 !This section has been created automatically by the script Abilint (TD).
2936 !Do not modify the following lines by hand.
2937 #undef ABI_FUNC
2938 #define ABI_FUNC 'getCoeffFromList'
2939 !End of the abilint section
2940 
2941  implicit none
2942 
2943 !Arguments ------------------------------------
2944 !scalar
2945  integer,intent(in) :: ia,ib,irpt,mu,ncoeff
2946  integer :: coeff
2947 !arrays
2948  integer,intent(in) :: list_coeff(6,ncoeff)
2949 !Local variables-------------------------------
2950 !scalar
2951  integer :: icoeff
2952 !arrays
2953 
2954 ! *************************************************************************
2955  coeff = 0
2956  do icoeff = 1,ncoeff
2957    if(mu==list_coeff(1,icoeff).and.&
2958 &     ia==list_coeff(2,icoeff).and.&
2959 &     ib==list_coeff(3,icoeff).and.&
2960 &     irpt==list_coeff(4,icoeff))then!.and.&
2961 !&     abs(weight-list_coeff(5,icoeff)) < tol16) then
2962      coeff = icoeff
2963      exit
2964    end if
2965  end do
2966 
2967 end function getCoeffFromList

m_polynomial_coeff/polynomial_coeff_broadcast [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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.

PARENTS

      m_effective_potential_file,m_fit_polynomial_coeff,m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

565 subroutine polynomial_coeff_broadcast(coefficients, source, comm)
566 
567 
568 !This section has been created automatically by the script Abilint (TD).
569 !Do not modify the following lines by hand.
570 #undef ABI_FUNC
571 #define ABI_FUNC 'polynomial_coeff_broadcast'
572 !End of the abilint section
573 
574  implicit none
575 
576 !Arguments ------------------------------------
577 !array
578  type(polynomial_coeff_type),intent(inout) :: coefficients
579  integer, intent(in) :: source,comm
580 
581 !Local variables-------------------------------
582 !scalars
583  integer :: ierr,ii
584 !arrays
585 
586 ! *************************************************************************
587 
588  if (xmpi_comm_size(comm) == 1) return
589 
590 ! Free the output
591  if (xmpi_comm_rank(comm) /= source) then
592    call polynomial_coeff_free(coefficients)
593  end if
594 
595  ! Transmit variables
596   call xmpi_bcast(coefficients%name, source, comm, ierr)
597   call xmpi_bcast(coefficients%nterm, source, comm, ierr)
598   call xmpi_bcast(coefficients%coefficient, source, comm, ierr)
599 
600  !Allocate arrays on the other nodes.
601   if (xmpi_comm_rank(comm) /= source) then
602     ABI_DATATYPE_ALLOCATE(coefficients%terms,(coefficients%nterm))
603     do ii=1,coefficients%nterm
604       call polynomial_term_free(coefficients%terms(ii))
605     end do
606   end if
607 ! Set the number of term on each node (needed for allocations of array)
608   do ii = 1,coefficients%nterm
609     call xmpi_bcast(coefficients%terms(ii)%ndisp, source, comm, ierr)
610     call xmpi_bcast(coefficients%terms(ii)%nstrain, source, comm, ierr)
611   end do
612 
613 ! Allocate arrays on the other nodes
614   if (xmpi_comm_rank(comm) /= source) then
615     do ii = 1,coefficients%nterm
616       ABI_ALLOCATE(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
617       coefficients%terms(ii)%atindx = 0
618       ABI_ALLOCATE(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
619       ABI_ALLOCATE(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
620       ABI_ALLOCATE(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
621       ABI_ALLOCATE(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
622       ABI_ALLOCATE(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
623     end do
624   end if
625 
626 ! Transfert value
627   do ii = 1,coefficients%nterm
628       call xmpi_bcast(coefficients%terms(ii)%weight, source, comm, ierr)
629       call xmpi_bcast(coefficients%terms(ii)%atindx, source, comm, ierr)
630       call xmpi_bcast(coefficients%terms(ii)%direction, source, comm, ierr)
631       call xmpi_bcast(coefficients%terms(ii)%cell, source, comm, ierr)
632       call xmpi_bcast(coefficients%terms(ii)%power_disp, source, comm, ierr)
633       call xmpi_bcast(coefficients%terms(ii)%power_strain, source, comm, ierr)
634       call xmpi_bcast(coefficients%terms(ii)%strain, source, comm, ierr)
635   end do
636 
637 
638 end subroutine polynomial_coeff_broadcast

m_polynomial_coeff/polynomial_coeff_evaluate [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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
  fcart(3,natom) = contribution  to the forces
  strten(6) = contribution to the stress tensor

PARENTS

      m_effective_potential

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

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

m_polynomial_coeff/polynomial_coeff_free [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_anharmonics_terms,m_effective_potential_file,m_fit_polynomial_coeff
      m_polynomial_coeff,mover_effpot

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

256 subroutine polynomial_coeff_free(polynomial_coeff)
257 
258 
259 !This section has been created automatically by the script Abilint (TD).
260 !Do not modify the following lines by hand.
261 #undef ABI_FUNC
262 #define ABI_FUNC 'polynomial_coeff_free'
263 !End of the abilint section
264 
265  implicit none
266 
267 !Arguments ------------------------------------
268 !scalars
269 !arrays
270  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
271 !Local variables-------------------------------
272 !scalar
273  integer :: ii
274 !arrays
275 
276 ! *************************************************************************
277 
278  if(allocated(polynomial_coeff%terms))then
279    do ii = 1,polynomial_coeff%nterm
280      call polynomial_term_free(polynomial_coeff%terms(ii))
281    end do
282    ABI_DATATYPE_DEALLOCATE(polynomial_coeff%terms)
283  end if
284  polynomial_coeff%name = ""
285  polynomial_coeff%nterm = 0
286  polynomial_coeff%coefficient = zero
287 
288 end subroutine polynomial_coeff_free

m_polynomial_coeff/polynomial_coeff_getList [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

1309 subroutine polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1310 &                                   natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,sc_size)
1311 
1312 
1313 !This section has been created automatically by the script Abilint (TD).
1314 !Do not modify the following lines by hand.
1315 #undef ABI_FUNC
1316 #define ABI_FUNC 'polynomial_coeff_getList'
1317 !End of the abilint section
1318 
1319  implicit none
1320 
1321 !Arguments ------------------------------------
1322 !scalars
1323  integer,intent(in) :: natom,nrpt
1324  integer,intent(out) :: ncoeff_sym,nstr_sym
1325 !arrays
1326  integer,intent(in) :: cell(3,nrpt)
1327  real(dp),intent(in):: dist(3,natom,natom,nrpt)
1328  type(crystal_t), intent(in) :: crystal
1329  integer,allocatable,intent(out) :: list_symcoeff(:,:,:),list_symstr(:,:,:)
1330  integer,optional,intent(in) :: sc_size(3)
1331  real(dp),intent(in):: range_ifc(3)
1332 !Local variables-------------------------------
1333 !scalar
1334  integer :: ia,ib,icoeff,icoeff2,icoeff_tot,icoeff_tmp,idisy1,idisy2,ii
1335  integer :: ipesy1,ipesy2,isym,irpt,irpt3,irpt_ref,irpt_sym
1336  integer :: jj,jsym,mu
1337  integer :: ncoeff,ncoeff2,ncoeff_max,nu
1338  integer :: nsym,shift_atm1(3)
1339  integer :: shift_atm2(3)
1340  real(dp):: dist_orig,dist_sym,tolsym8
1341  logical :: found,check_pbc,possible
1342 !arrays
1343  integer :: isym_rec(3,3),isym_rel(3,3),sc_size_in(3)
1344  integer :: transl(3),min_range(3),max_range(3)
1345  integer,allocatable :: blkval(:,:,:,:,:),list(:),list_symcoeff_tmp(:,:,:),list_symcoeff_tmp2(:,:,:)
1346  integer,allocatable :: list_symstr_tmp(:,:,:),indsym(:,:,:) ,symrec(:,:,:),symrel(:,:,:)
1347  real(dp),allocatable :: tnons(:,:)
1348  real(dp),allocatable :: wkdist(:),xcart(:,:),xred(:,:),distance(:,:,:)
1349  real(dp) :: difmin(3)
1350  real(dp) :: rprimd(3,3)
1351  real(dp) :: tratom(3)
1352  character(len=500) :: message
1353 
1354 ! *************************************************************************
1355 
1356 
1357 !Initialisation of variables
1358  irpt_sym = 0
1359  nsym   = crystal%nsym
1360  rprimd = crystal%rprimd
1361  ABI_ALLOCATE(xcart,(3,natom))
1362  ABI_ALLOCATE(xred,(3,natom))
1363  xcart(:,:) = crystal%xcart(:,:)
1364  xred(:,:)  = crystal%xred(:,:)
1365  ncoeff_max = nrpt*natom*natom*3*3
1366 
1367 !Found the ref cell
1368  irpt_ref = 1
1369  do irpt=1,nrpt
1370    if(all(cell(:,irpt)==0))then
1371      irpt_ref = irpt
1372 !     exit
1373    end if
1374  end do
1375 
1376  !Set the size of the interaction
1377  check_pbc = .FALSE.
1378  sc_size_in = 0
1379  min_range = 0; max_range = 0
1380  if(present(sc_size))then
1381    sc_size_in = sc_size
1382    do mu=1,3
1383      call findBound_supercell(min_range(mu),max_range(mu),sc_size_in(mu))
1384    end do
1385  end if
1386 
1387 !Obtain a list of rotated atom labels:
1388  ABI_ALLOCATE(indsym,(4,nsym,natom))
1389  ABI_ALLOCATE(symrec,(3,3,nsym))
1390  ABI_ALLOCATE(symrel,(3,3,nsym))
1391  ABI_ALLOCATE(tnons,(3,nsym))
1392  symrec = crystal%symrec
1393  symrel = crystal%symrel
1394  tnons  = crystal%tnons
1395 
1396  tolsym8=tol20
1397  call symatm(indsym,natom,nsym,symrec,tnons,&
1398 &            tolsym8,crystal%typat,crystal%xred)
1399  ABI_ALLOCATE(blkval,(3,natom,3,natom,nrpt))
1400  ABI_ALLOCATE(list,(natom*nrpt))
1401  ABI_ALLOCATE(list_symcoeff_tmp,(5,ncoeff_max,nsym))
1402  ABI_ALLOCATE(wkdist,(natom*nrpt))
1403 
1404 !1-Fill strain list
1405  ABI_ALLOCATE(list_symstr_tmp,(6,nsym,2))
1406  list_symstr_tmp = 1
1407  do ia=1,6
1408    if(list_symstr_tmp(ia,1,1)==0)cycle
1409 !  Transform the voigt notation
1410    if(ia<=3)then
1411      mu=ia;nu=ia
1412    else
1413      select case(ia)
1414      case(4)
1415        mu=2;nu=3
1416      case(5)
1417        mu=1;nu=3
1418      case(6)
1419        mu=1;nu=2
1420      end select
1421    end if
1422    do isym=1,nsym
1423 !    Get the symmetry matrix
1424      isym_rel(:,:) = crystal%symrel(:,:,isym)
1425      do idisy1=1,3
1426        do idisy2=1,3
1427          if((isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)) then
1428 !          Transform to the voig notation
1429            if(idisy1==idisy2)then
1430              list_symstr_tmp(ia,isym,1) = idisy1
1431              list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1)
1432            else
1433              if(idisy1==1.or.idisy2==1)then
1434                if(idisy1==2.or.idisy2==2)then
1435                  list_symstr_tmp(ia,isym,1) = 6
1436                end if
1437                if(idisy1==3.or.idisy2==3)then
1438                  list_symstr_tmp(ia,isym,1) = 5
1439                end if
1440              else
1441                list_symstr_tmp(ia,isym,1) = 4
1442              end if
1443            end if
1444            list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1) * isym_rel(nu,idisy2)
1445          end if
1446        end do
1447      end do
1448 !    Remove the symetric
1449 !     if(list_symstr_tmp(ia,isym,1) > ia) then
1450 !       list_symstr_tmp(list_symstr_tmp(ia,isym,1),:,1) = 0
1451 !     end if
1452    end do
1453  end do
1454 
1455 !Count the number of strain and transfert into the final array
1456   nstr_sym = 0
1457   do ia=1,6
1458     if(list_symstr_tmp(ia,1,1)/=0) nstr_sym = nstr_sym + 1
1459   end do
1460 
1461  if(allocated(list_symstr))then
1462    ABI_DEALLOCATE(list_symstr)
1463  end if
1464  ABI_ALLOCATE(list_symstr,(nstr_sym,nsym,2))
1465 
1466  icoeff_tmp = 1
1467  do ia=1,6
1468    if(list_symstr_tmp(ia,1,1)/=0) then
1469      list_symstr(icoeff_tmp,:,:) = list_symstr_tmp(ia,:,:)
1470      icoeff_tmp = icoeff_tmp + 1
1471    end if
1472  end do
1473 !END STRAIN
1474 
1475 !Compute the distance between each atoms. Indeed the dist array contains the difference of
1476 !cartesian coordinate for each direction
1477  ABI_ALLOCATE(distance,(natom,natom,nrpt))
1478  do ia=1,natom
1479    do ib=1,natom
1480      do irpt=1,nrpt
1481        distance(ia,ib,irpt) = ((dist(1,ia,ib,irpt))**2+(dist(2,ia,ib,irpt))**2+&
1482 &                              (dist(3,ia,ib,irpt))**2)**0.5
1483      end do
1484    end do
1485  end do
1486 
1487 
1488 !Set to one blkval, all the coeff have to be compute
1489  blkval = 1
1490  icoeff = 1
1491  icoeff_tot = 1
1492  list_symcoeff_tmp = 0
1493 
1494 !2-Fill atom list
1495 !Big loop over generic atom
1496  do ia=1,natom
1497    wkdist(:)=reshape(distance(ia,:,:),(/natom*nrpt/))
1498    do ii=1,natom*nrpt
1499      list(ii)=ii
1500    end do
1501    call sort_dp(natom*nrpt,wkdist,list,tol8)
1502    do ii=1,natom*nrpt
1503 !    Get the irpt and ib
1504      irpt=(list(ii)-1)/natom+1
1505      ib=list(ii)-natom*(irpt-1)
1506      possible = .true.
1507 !Old way with the cut off
1508 !     if(((dist(1,ia,ib,irpt)**2+dist(2,ia,ib,irpt)**2+dist(3,ia,ib,irpt)**2)**0.5) > 9)then
1509 !       possible = .false.
1510 !     end if
1511      do jj=1,3
1512 !        if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)  > tol10.or.&
1513 ! &          abs(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj))  < tol10)then
1514         if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)  > tol10)then
1515        possible = .false.
1516        end if
1517      end do
1518 
1519 !    If this distance is superior to the cutoff, we don't compute that term
1520      if(.not.possible)then
1521        blkval(:,ia,:,ib,irpt)= 0
1522        if(irpt==irpt_ref)blkval(:,ib,:,ia,irpt)= 0
1523 !        Stop the loop
1524        cycle
1525      end if
1526 
1527 !    If this coefficient is not possible, we cycle...
1528      if (all(blkval(:,ia,:,ib,irpt)==0)) cycle
1529 
1530 !    Save the distance between the two atoms for futur checks
1531      dist_orig = (dist(1,ia,ib,irpt)**2+dist(2,ia,ib,irpt)**2+dist(3,ia,ib,irpt)**2)**0.5
1532 
1533      do mu=1,3
1534        do nu=1,3
1535 !      Check if : - The coefficient is not yet compute
1536 !                 - The directions are the same
1537 !                 - The atoms are not equivalent
1538          if (mu/=nu) then
1539            blkval(mu,ia,nu,ib,irpt)=0
1540            blkval(nu,ia,mu,ib,irpt)=0
1541            cycle
1542          end if
1543 !        Pass if the atoms are identical and in the ref cell
1544          if(irpt==irpt_ref.and.ia==ib) then
1545            blkval(mu,ia,nu,ib,irpt)=0
1546            blkval(nu,ib,mu,ia,irpt)=0
1547            cycle
1548          end if
1549          if(blkval(mu,ia,nu,ib,irpt)==1)then
1550 !          Loop over symmetries
1551            do isym=1,nsym
1552 !            Get the symmetry matrix for this sym
1553              isym_rec(:,:)  = crystal%symrec(:,:,isym)
1554              isym_rel(:,:) = crystal%symrel(:,:,isym)
1555 !            Get the corresponding atom and shift with the symetries
1556 !            For atom 1
1557              ipesy1 = indsym(4,isym,ia)
1558              shift_atm1 = indsym(1:3,isym,ia)
1559 !            And atom 2
1560              do jj=1,3 ! Apply transformation to original coordinates.
1561               tratom(jj) = dble(isym_rec(1,jj))*(xred(1,ib)+cell(1,irpt)-tnons(1,isym))&
1562 &                         +dble(isym_rec(2,jj))*(xred(2,ib)+cell(2,irpt)-tnons(2,isym))&
1563 &                         +dble(isym_rec(3,jj))*(xred(3,ib)+cell(3,irpt)-tnons(3,isym))
1564 
1565              end do
1566 
1567 !            Find symmetrically equivalent atom
1568              call symchk(difmin,ipesy2,natom,tratom,transl,crystal%typat(ib),&
1569 &                        crystal%typat,xred(:,:))
1570 
1571 !            Put information into array indsym: translations and label
1572              shift_atm2(:)= transl(:) - shift_atm1(:)
1573              found = .false.
1574              do irpt3=1,nrpt
1575                if(cell(1,irpt3)==shift_atm2(1).and.&
1576 &                 cell(2,irpt3)==shift_atm2(2).and.&
1577 &                 cell(3,irpt3)==shift_atm2(3))then
1578                  found = .true.
1579                  irpt_sym = irpt3
1580                end if
1581              end do
1582 
1583 !            Check the distance
1584              dist_sym = (dist(1,ipesy1,ipesy2,irpt_sym)**2+&
1585 &                        dist(2,ipesy1,ipesy2,irpt_sym)**2+&
1586 &                        dist(3,ipesy1,ipesy2,irpt_sym)**2)**0.5
1587              if(abs(dist_orig - dist_sym) > tol10)then
1588                write(message, '(a,i0,2a,I0,a,es15.8,2a,es15.8,2a)' )&
1589 &                'The distance between the atoms for the coefficient number ',icoeff,ch10,&
1590 &                'with the symmetry ',isym,' is ',dist_sym,ch10,'but the original distance is',&
1591 &                   dist_orig,ch10,&
1592 &                'Action: Contact abinit group'
1593                MSG_BUG(message)
1594              end if
1595 !            Now that a symmetric perturbation has been obtained,
1596 !            including the expression of the symmetry matrix, see
1597 !            if the symmetric perturbations are available
1598              do idisy1=1,3
1599                do idisy2=1,3
1600                  if (idisy1/=idisy2) then
1601 !                  Remove this term (is not computed)
1602 !                  Also remove opposite term... (Srx-Tix) = (Ti-Srx)
1603                    blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1604                    blkval(idisy2,ipesy1,idisy1,ipesy2,irpt_sym) = 0
1605                    cycle
1606                  else
1607                    if(isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)then
1608                      if(.not.found.or.(irpt_sym==irpt_ref.and.ipesy1==ipesy2)) then
1609 !                      Remove this term (is not computed) Sr-Sr or not include in the cell
1610 !                      Also remove oposite term... (Srx-Tix) = (Ti-Srx)
1611                        blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0
1612                        blkval(idisy2,ipesy2,idisy1,ipesy1,irpt_sym) = 0
1613                        cycle
1614                      else
1615 !                      Fill the list with the coeff and symmetric (need all symmetrics)
1616                        list_symcoeff_tmp(1:4,icoeff,isym)=(/idisy1,ipesy1,ipesy2,irpt_sym/)
1617 !                      Check the sign
1618                        if(isym_rel(mu,idisy1)/=isym_rel(nu,idisy2))then
1619                          write(message, '(a,i0,a,I0,4a)' )&
1620 &                        'The sign of coefficient number ',icoeff,' with the symmetry ',isym,ch10,&
1621 &                        'can not be found... Something is going wrong',ch10,&
1622 &                        'Action: Contact abinit group'
1623                          MSG_BUG(message)
1624                        end if
1625                        list_symcoeff_tmp(5,icoeff,isym)= isym_rel(nu,idisy2)
1626                      end if
1627                    end if
1628                  end if
1629                end do
1630              end do
1631            end do ! end loop sym
1632            icoeff = icoeff + 1
1633          end if
1634 !        This coeff is now computed
1635          blkval(mu,ia,nu,ib,irpt)= 0
1636        end do ! end loop nu
1637      end do ! end loop mu
1638    end do ! end loop ii
1639  end do ! end loop ia
1640 
1641 !Reset the output
1642  if(allocated(list_symcoeff))then
1643    ABI_DEALLOCATE(list_symcoeff)
1644  end if
1645 
1646  ABI_DEALLOCATE(distance)
1647 
1648 !Transfert the final array with all the coefficients
1649 !With this array, we can access to all the terms presents
1650 !ncoeff1 + symetrics
1651 !first dimension is 4 (mu,ia,ib,irpt)
1652 !irpt is the index of the cell of the atom ib in the cell array
1653 !example cell(:,irpt=12) can be (-1 0 -2). The cell of ia is
1654 !always 0 0 0
1655 !Transfert the final array for the list of irreductible coeff and symetries
1656 !With this array, we can access to the irretuctible coefficients (ncoeff1) and
1657 !all the symetrics of these coefficients (nsym)
1658 !first dimension is 5 (mu,ia,ib,irpt,icoeff)
1659 !icoeff is the position of this coefficients in the list_fullcoeff array
1660 
1661 !1/ step remove the zero coeff in this array
1662  ncoeff = 0
1663  do icoeff = 1,ncoeff_max
1664      if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then
1665      ncoeff = ncoeff + 1
1666    end if
1667  end do
1668 
1669  ABI_ALLOCATE(list_symcoeff_tmp2,(6,ncoeff,nsym))
1670  list_symcoeff_tmp2 = 0
1671  icoeff = 0
1672  do icoeff_tmp = 1,ncoeff_max
1673    if(.not.(all(list_symcoeff_tmp(:,icoeff_tmp,1)==0)))then
1674      icoeff = icoeff + 1
1675      list_symcoeff_tmp2(1:5,icoeff,:) = list_symcoeff_tmp(1:5,icoeff_tmp,:)
1676      end if
1677  end do
1678 
1679 
1680 !2/ set the dimension six of list_symcoeff_tmp2(6,icoeffs,1)
1681 !   and check is a symetric coeff is not coresspondig to an other
1682 !   one, in this case we set this coeff to 0
1683 ! ncoeff2 = zero
1684  do icoeff = 1,ncoeff
1685 !  found the index of each coeff in list_fullcoeff
1686    do isym = 1,nsym
1687      icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,1),&
1688 &                               list_symcoeff_tmp2(2,icoeff,isym),&
1689 &                               list_symcoeff_tmp2(3,icoeff,isym),&
1690 &                               list_symcoeff_tmp2(4,icoeff,isym),&
1691 &                               list_symcoeff_tmp2(1,icoeff,isym),&
1692 &                               ncoeff)
1693      list_symcoeff_tmp2(6,icoeff,isym) = icoeff2
1694    end do
1695  end do
1696 
1697 !2.5/do checks
1698  do icoeff = 1,ncoeff
1699    do isym = 1,nsym
1700      if(list_symcoeff_tmp2(6,icoeff,isym)==0)then
1701        write(message, '(a,i0,a,I0,4a)' )&
1702 &           'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1703 &           'has no equivalent',ch10,&
1704 &           'Action: Contact abinit group'
1705        MSG_BUG(message)
1706      else
1707        if(icoeff /= list_symcoeff_tmp2(6,icoeff,isym))then
1708          if(list_symcoeff_tmp2(1,icoeff,isym)/=&
1709 &           list_symcoeff_tmp2(1,list_symcoeff_tmp2(6,icoeff,isym),1))then
1710            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1711 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1712 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1713 &          'because the direction is different:',ch10,&
1714 &          'Action: Contact abinit group'
1715            MSG_BUG(message)
1716          end if
1717          if(list_symcoeff_tmp2(4,icoeff,isym)/=&
1718 &           list_symcoeff_tmp2(4,list_symcoeff_tmp2(6,icoeff,isym),1))then
1719            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1720 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1721 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1722 &          'because the cell is different',ch10,&
1723 &          'Action: Contact abinit group'
1724            MSG_BUG(message)
1725          end if
1726          if((list_symcoeff_tmp2(2,icoeff,isym)/=&
1727 &            list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1728 &            list_symcoeff_tmp2(3,icoeff,isym)/=&
1729 &            list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1)).and.&
1730 &           (list_symcoeff_tmp2(2,icoeff,isym)/=&
1731 &            list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1).and.&
1732 &            list_symcoeff_tmp2(3,icoeff,isym)/=&
1733 &            list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1)))then
1734            write(message, '(a,i0,a,I0,2a,I0,4a)' )&
1735 &          'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,&
1736 &          'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,&
1737 &          'because the atoms different',ch10,&
1738 &          'Action: Contact abinit group'
1739            MSG_BUG(message)
1740          end if
1741        end if
1742      end if
1743    end do
1744  end do
1745 
1746 
1747 !Check if the atom 2 is not (with the PBC) is in the same cell than
1748 !the atom 1. For example if you want to fit 2x2x2 cell the interation
1749 !Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process...
1750 !If check_pbc is true we remove these kind of terms
1751  do icoeff = 1,ncoeff
1752    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1753      do mu=1,3
1754        if(min_range(mu) > cell(mu,list_symcoeff_tmp2(4,icoeff,1)) .or. &
1755 &       cell(mu,list_symcoeff_tmp2(4,icoeff,1)) > max_range(mu))then
1756          list_symcoeff_tmp2(:,icoeff,:)=0
1757          exit
1758        end if
1759      end do
1760    end if
1761  end do
1762 
1763 !3/ Remove useless terms like opposites
1764  do icoeff = 1,ncoeff
1765    do isym = 1,nsym
1766      icoeff2 = list_symcoeff_tmp2(6,icoeff,isym)
1767      if (icoeff2> icoeff)then
1768        list_symcoeff_tmp2(:,icoeff2,1) = 0
1769      end if
1770      do jsym=1,nsym
1771        icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),&
1772 &                                 list_symcoeff_tmp2(3,icoeff,isym),&
1773 &                                 list_symcoeff_tmp2(2,icoeff,isym),&
1774 &                                 list_symcoeff_tmp2(4,icoeff,isym),&
1775 &                                 list_symcoeff_tmp2(1,icoeff,isym),&
1776 &                                 ncoeff)
1777        if (icoeff2> icoeff)then
1778          list_symcoeff_tmp2(:,icoeff2,1) = 0
1779        end if
1780      end do
1781    end do
1782  end do
1783 
1784 !4/ Recount the number of coeff after step 3
1785  ncoeff2 = 0
1786  do icoeff = 1,ncoeff
1787    if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then
1788      ncoeff2 = ncoeff2 + 1
1789    end if
1790  end do
1791 
1792  ABI_DEALLOCATE(list_symcoeff_tmp)
1793  ABI_ALLOCATE(list_symcoeff_tmp,(6,ncoeff2,nsym))
1794 
1795  list_symcoeff_tmp = 0
1796  icoeff = 0
1797 
1798  do icoeff_tmp = 1,ncoeff
1799    if(.not.(all(list_symcoeff_tmp2(:,icoeff_tmp,1)==0)))then
1800      icoeff = icoeff + 1
1801      list_symcoeff_tmp(:,icoeff,:) = list_symcoeff_tmp2(:,icoeff_tmp,:)
1802    end if
1803  end do
1804 
1805 !5/ Final transfert
1806  ABI_ALLOCATE(list_symcoeff,(6,ncoeff2,nsym))
1807  list_symcoeff = 0
1808  icoeff = 0
1809  do icoeff = 1,ncoeff2
1810    list_symcoeff(1:6,icoeff,:) = list_symcoeff_tmp(1:6,icoeff,:)
1811    do isym=1,nsym
1812    end do
1813  end do
1814 
1815 !6/ reset 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  do icoeff = 1,ncoeff2
1819 !  found the index of each coeff in list_fullcoeff
1820    do isym = 1,nsym
1821      icoeff2 = getCoeffFromList(list_symcoeff(:,:,1),&
1822 &                               list_symcoeff(2,icoeff,isym),&
1823 &                               list_symcoeff(3,icoeff,isym),&
1824 &                               list_symcoeff(4,icoeff,isym),&
1825 &                               list_symcoeff(1,icoeff,isym),&
1826 &                               ncoeff)
1827      list_symcoeff(6,icoeff,isym) = icoeff2
1828    end do
1829  end do
1830 
1831 !Set the max number of coeff inside list_symcoeff
1832  ncoeff_sym = ncoeff2
1833 
1834 !Deallocation
1835  ABI_DEALLOCATE(blkval)
1836  ABI_DEALLOCATE(list)
1837  ABI_DEALLOCATE(list_symcoeff_tmp)
1838  ABI_DEALLOCATE(list_symcoeff_tmp2)
1839  ABI_DEALLOCATE(list_symstr_tmp)
1840  ABI_DEALLOCATE(indsym)
1841  ABI_DEALLOCATE(symrec)
1842  ABI_DEALLOCATE(symrel)
1843  ABI_DEALLOCATE(tnons)
1844  ABI_DEALLOCATE(xcart)
1845  ABI_DEALLOCATE(xred )
1846  ABI_DEALLOCATE(wkdist)
1847 
1848 end subroutine polynomial_coeff_getList

m_polynomial_coeff/polynomial_coeff_getName [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_effective_potential_file,m_fit_polynomial_coeff,m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

417 subroutine polynomial_coeff_getName(name,polynomial_coeff,symbols,recompute,iterm)
418 
419 
420 !This section has been created automatically by the script Abilint (TD).
421 !Do not modify the following lines by hand.
422 #undef ABI_FUNC
423 #define ABI_FUNC 'polynomial_coeff_getName'
424 !End of the abilint section
425 
426  implicit none
427 
428 !Arguments ------------------------------------
429 !scalars
430  integer,optional,intent(in) :: iterm
431 !arrays
432  character(len=5),intent(in) :: symbols(:)
433  character(len=200),intent(out):: name
434  type(polynomial_coeff_type),optional, intent(in) :: polynomial_coeff
435  logical,optional,intent(in) :: recompute
436 !Local variables-------------------------------
437 !scalar
438  integer :: ii,idisp,iterm_in
439  logical :: need_recompute = .FALSE.
440 !arrays
441  integer :: cell_atm1(3),cell_atm2(3)
442  character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
443  character(len=1) :: dir
444  character(len=2) :: power_disp,power_dispchar
445  character(len=5) :: atm1,atm2
446  character(len=100):: atm1_tmp,atm2_tmp
447  character(len=200):: text
448  character(len=500):: message
449 ! *************************************************************************
450 
451 !Reset output
452  name=""
453  iterm_in = 1
454 
455 !Set the optional arguments
456  if(present(recompute)) need_recompute = recompute
457  if(present(iterm)) then
458    iterm_in = iterm
459  else
460    if(need_recompute)then
461      iterm_in = -1
462      do ii=1,polynomial_coeff%nterm
463 !      Find the index of the ref
464        if(iterm_in==-1) then !Need to find the reference term
465          do idisp=1,polynomial_coeff%terms(ii)%ndisp
466            if(polynomial_coeff%terms(ii)%direction(idisp) > 0) then
467              iterm_in = ii
468              if(any(polynomial_coeff%terms(ii)%cell(:,1,idisp) /= 0).or.&
469 &               any(polynomial_coeff%terms(ii)%cell(:,2,idisp) /= 0)) then
470                iterm_in = -1
471                exit
472              end if
473            end if
474          end do!end do disp
475        else
476          exit
477        end if
478      end do!end do term
479 !    If not find, we set to the first element
480      if(iterm_in==-1) iterm_in = 1
481    else
482      iterm_in = 1
483    end if
484  end if
485 
486 !Do check
487  if(iterm_in > polynomial_coeff%nterm.or.iterm_in < 0) then
488    write(message, '(5a)')&
489 &      ' The number of the requested term for the generation of',ch10,&
490 &      'the name of the coefficient is not possible.',ch10,&
491 &      'Action: Contact Abinit group.'
492    MSG_BUG(message)
493  end if
494 
495  if(polynomial_coeff%name /= "".and..not.need_recompute)then
496    name = polynomial_coeff%name
497  else
498 !  Nedd to recompute
499    do idisp=1,polynomial_coeff%terms(iterm_in)%ndisp
500      text = ""
501      !Fill variables for this displacement
502      write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_disp(idisp)
503      power_disp=trim(power_dispchar)
504 
505      atm1=symbols(polynomial_coeff%terms(iterm_in)%atindx(1,idisp))
506      atm2=symbols(polynomial_coeff%terms(iterm_in)%atindx(2,idisp))
507      dir=mutodir(polynomial_coeff%terms(iterm_in)%direction(idisp))
508      cell_atm1=polynomial_coeff%terms(iterm_in)%cell(:,1,idisp)
509      cell_atm2=polynomial_coeff%terms(iterm_in)%cell(:,2,idisp)
510 !    Construct ATM1
511      if (any(cell_atm1(:) /= 0) )then
512        write(atm1_tmp,'(4a,I0,a,I0,a,I0,a)')  trim(atm1),"_",dir,"[",cell_atm1(1)," ",&
513 &                                               cell_atm1(2)," ",cell_atm1(3),"]"
514      else
515        atm1_tmp = trim(atm1)//"_"//dir
516      end if
517 !      Construct ATM2
518      if(any(cell_atm2(:) /= 0))then
519        write(atm2_tmp,'(4a,I0,a,I0,a,I0,a)')  trim(atm2),"_",dir,"[",cell_atm2(1)," ",&
520  &                                              cell_atm2(2)," ",cell_atm2(3),"]"
521      else
522        atm2_tmp = trim(atm2)//"_"//dir
523      end if
524 
525      text="("//trim(atm1_tmp)//"-"//trim(atm2_tmp)//")^"//power_disp
526      name = trim(name)//trim(text)
527    end do
528    !Strain case
529    do idisp=1,polynomial_coeff%terms(iterm_in)%nstrain
530      write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_strain(idisp)
531      power_disp=trim(power_dispchar)
532      dir=mutodir(3+polynomial_coeff%terms(iterm_in)%strain(idisp))
533      text="("//"eta_"//trim(dir)//")^"//power_disp
534      name = trim(name)//trim(text)
535    end do
536  end if
537 
538 end subroutine polynomial_coeff_getName

m_polynomial_coeff/polynomial_coeff_getNorder [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_fit_polynomial_coeff,mover_effpot

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

1895 subroutine polynomial_coeff_getNorder(coefficients,crystal,cutoff,ncoeff,ncoeff_tot,power_disps,&
1896 &                                     max_power_strain,option,sc_size,comm,anharmstr,spcoupling,&
1897 &                                     distributed,only_odd_power,only_even_power,verbose)
1898 
1899 
1900 !This section has been created automatically by the script Abilint (TD).
1901 !Do not modify the following lines by hand.
1902 #undef ABI_FUNC
1903 #define ABI_FUNC 'polynomial_coeff_getNorder'
1904 !End of the abilint section
1905 
1906  implicit none
1907 
1908 !Arguments ------------------------------------
1909 !scalars
1910  integer,intent(in) :: max_power_strain,option,comm
1911  integer,intent(out):: ncoeff,ncoeff_tot
1912  real(dp),intent(in):: cutoff
1913  logical,optional,intent(in) :: anharmstr,spcoupling,distributed,verbose
1914  logical,optional,intent(in) :: only_odd_power,only_even_power
1915 !arrays
1916  integer,intent(in) :: power_disps(2),sc_size(3)
1917  type(crystal_t), intent(inout) :: crystal
1918  type(polynomial_coeff_type),allocatable,intent(inout) :: coefficients(:)
1919 !Local variables-------------------------------
1920 !scalar
1921  integer :: ia,ib,icoeff,icoeff2,icoeff3,ierr,ii,irpt,irpt_ref,iterm
1922  integer :: lim1,lim2,lim3
1923  integer :: master,my_rank,my_ncoeff,my_newncoeff,natom,ncombination,ncoeff_max,ncoeff_sym
1924  integer :: ncoeff_alone,ndisp_max,nproc,nrpt,nsym,nterm,nstr_sym,r1,r2,r3,my_size
1925  integer :: my_icoeff,rank_to_send,rank_to_receive,rank_to_send_save
1926  real(dp):: norm
1927  logical :: iam_master,need_anharmstr,need_spcoupling,need_distributed,need_verbose
1928  logical :: need_only_odd_power,need_only_even_power
1929 !arrays
1930  integer :: ncell(3)
1931  integer,allocatable :: buffsize(:),buffdispl(:)
1932  integer,allocatable :: cell(:,:),compatibleCoeffs(:,:)
1933  integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:),list_coeff(:),list_combination(:,:)
1934  integer,allocatable :: list_combination_tmp(:,:)
1935  integer,allocatable  :: my_coefflist(:),my_coeffindexes(:),my_newcoeffindexes(:)
1936  real(dp) :: rprimd(3,3),range_ifc(3)
1937  real(dp),allocatable :: dist(:,:,:,:),rpt(:,:)
1938  real(dp),allocatable :: xcart(:,:),xred(:,:)
1939  character(len=5),allocatable :: symbols(:)
1940  character(len=200):: name
1941  character(len=500) :: message
1942  type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_tmp
1943  type(polynomial_term_type),dimension(:),allocatable :: terms
1944 ! character(len=fnlen) :: filename
1945 ! *************************************************************************
1946 
1947  !MPI variables
1948  master = 0
1949  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1950  iam_master = (my_rank == master)
1951 
1952 !Free the output
1953  if(allocated(coefficients))then
1954    do ii =1,size(coefficients)
1955      call polynomial_coeff_free(coefficients(ii))
1956    end do
1957    ABI_DEALLOCATE(coefficients)
1958  end if
1959 
1960 !Check
1961  if(option > power_disps(2))then
1962    write(message, '(3a)' )&
1963 &       'Option can not be superior to the maximum order ',ch10,&
1964 &       'Action: contact abinit group'
1965    MSG_ERROR(message)
1966  end if
1967 
1968 !Initialisation of variables
1969  need_anharmstr = .TRUE.
1970  if(present(anharmstr)) need_anharmstr = anharmstr
1971  need_spcoupling = .TRUE.
1972  if(present(spcoupling)) need_spcoupling = spcoupling
1973  need_distributed = .FALSE.
1974  if(present(distributed)) need_distributed  = distributed
1975  need_verbose = .TRUE.
1976  if(present(verbose)) need_verbose = verbose
1977  need_only_odd_power = .FALSE.
1978  if(present(only_odd_power)) need_only_odd_power = only_odd_power
1979  need_only_even_power = .FALSE.
1980  if(present(only_even_power)) need_only_even_power = only_even_power
1981 
1982  if(need_only_odd_power.and.need_only_even_power)then
1983       write(message, '(3a)' )&
1984 &       'need_only_odd_power and need_only_even_power are both true',ch10,&
1985 &       'Action: contact abinit group'
1986    MSG_ERROR(message)
1987  end if
1988 
1989  natom  = crystal%natom
1990  nsym   = crystal%nsym
1991  rprimd = crystal%rprimd
1992 
1993  ABI_ALLOCATE(xcart,(3,natom))
1994  ABI_ALLOCATE(xred,(3,natom))
1995  xcart(:,:) = crystal%xcart(:,:)
1996  xred(:,:)  = crystal%xred(:,:)
1997 
1998 !Compute the max range of the ifc with respect to the trainning set
1999  range_ifc(:) = zero
2000  do ii=1,3
2001    norm = sqrt(rprimd(ii,1)**2+ rprimd(ii,2)**2+rprimd(ii,3)**2)
2002    range_ifc(ii) = range_ifc(ii) + norm * sc_size(ii) / 2.0
2003  end do
2004 
2005 !compute new ncell
2006  ncell = sc_size
2007  lim1=((ncell(1)/2)) + 1
2008  lim2=((ncell(2)/2)) + 1
2009  lim3=((ncell(3)/2)) + 1
2010  if(mod(ncell(1),2)/=0) lim1=lim1+1
2011  if(mod(ncell(2),2)/=0) lim2=lim2+1
2012  if(mod(ncell(3),2)/=0) lim3=lim3+1
2013  nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)
2014 
2015  ncell(1) = 2*lim1+1
2016  ncell(2) = 2*lim2+1
2017  ncell(3) = 2*lim3+1
2018 
2019  !Build the rpt point
2020  ABI_ALLOCATE(rpt,(3,nrpt))
2021  ABI_ALLOCATE(cell,(3,nrpt))
2022 
2023 !WARNING:
2024 !Put the reference cell into the first element
2025 !the code will first deal with the atoms of the first cell
2026  irpt = 1
2027  irpt_ref = 1
2028  rpt(:,1) = zero
2029  cell(:,irpt)=0
2030 !Fill other rpt:
2031  do r1=lim1,-lim1,-1
2032    do r2=lim2,-lim2,-1
2033      do r3=lim3,-lim3,-1
2034        if(r1==0.and.r2==0.and.r3==0) then
2035          cycle
2036        end if
2037        irpt=irpt+1
2038        rpt(1,irpt)=r1*rprimd(1,1)+r2*rprimd(1,2)+r3*rprimd(1,3)
2039        rpt(2,irpt)=r1*rprimd(2,1)+r2*rprimd(2,2)+r3*rprimd(2,3)
2040        rpt(3,irpt)=r1*rprimd(3,1)+r2*rprimd(3,2)+r3*rprimd(3,3)
2041        cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3
2042      end do
2043    end do
2044  end do
2045 
2046  ABI_ALLOCATE(symbols,(natom))
2047  call symbols_crystal(crystal%natom,crystal%ntypat,crystal%npsp,&
2048 &                     symbols,crystal%typat,crystal%znucl)
2049 
2050 !Compute the distances between atoms
2051 !Now dist(3,ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt.
2052  ABI_ALLOCATE(dist,(3,natom,natom,nrpt))
2053  dist = zero
2054  do ia=1,natom
2055    do ib=1,natom
2056      do irpt=1,nrpt
2057        dist(1,ia,ib,irpt) = xcart(1,ib)-xcart(1,ia)+rpt(1,irpt)
2058        dist(2,ia,ib,irpt) = xcart(2,ib)-xcart(2,ia)+rpt(2,irpt)
2059        dist(3,ia,ib,irpt) = xcart(3,ib)-xcart(3,ia)+rpt(3,irpt)
2060      end do
2061    end do
2062  end do
2063 
2064  if(need_verbose)then
2065    write(message,'(1a)')' Generation of the list of all the possible coefficients'
2066    call wrtout(std_out,message,'COLL')
2067  end if
2068  call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
2069 &                              natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,sc_size=sc_size)
2070 
2071  ABI_DEALLOCATE(dist)
2072  ABI_DEALLOCATE(rpt)
2073 
2074 !Compute the total number of coefficient
2075  ncoeff_tot = ncoeff_sym+nstr_sym
2076 
2077  if(iam_master)then
2078 !Check the distanceance bewteen coefficients and store integer:
2079 ! 0: the mix between these coefficient is not possible
2080 ! 1: the mix between these coefficient is possible
2081    ABI_ALLOCATE(compatibleCoeffs,(ncoeff_tot,ncoeff_tot))
2082    compatibleCoeffs(:,:) = 1
2083 
2084    if(need_verbose)then
2085      write(message,'(1a)')' Check the compatible coefficients with respect to the cutoff'
2086      call wrtout(std_out,message,'COLL')
2087    end if
2088 
2089    do icoeff=1,ncoeff_tot
2090      do icoeff2=1,ncoeff_tot
2091 !      Select case:
2092 !      if both icoeff are displacement => check the distance
2093 !      if both icoeff are strain => check the flag
2094 !      Otherwise cycle (we keep the term)
2095        if(icoeff>ncoeff_sym.and.icoeff2<=ncoeff_sym)cycle
2096        if(icoeff2<=ncoeff_sym.and.icoeff2>ncoeff_sym)cycle
2097        if((icoeff>ncoeff_sym.or.icoeff2>ncoeff_sym).and.&
2098 &       .not.need_anharmstr.and..not.need_spcoupling) then
2099          compatibleCoeffs(icoeff,icoeff2) = 0
2100        end if
2101        if(icoeff2<=ncoeff_sym.and.icoeff2<=ncoeff_sym)then
2102          if(distance_supercell(xcart(:,list_symcoeff(2,icoeff,1)),&
2103 &                  xcart(:,list_symcoeff(2,icoeff2,1)),rprimd,&
2104 &                  cell(:,1),cell(:,1))>=cutoff.or.&
2105 &         distance_supercell(xcart(:,list_symcoeff(2,icoeff,1)),&
2106 &                  xcart(:,list_symcoeff(3,icoeff2,1)),rprimd,&
2107 &                  cell(:,1),cell(:,list_symcoeff(4,icoeff2,1)))>=cutoff.or.&
2108 &         distance_supercell(xcart(:,list_symcoeff(3,icoeff,1)),&
2109 &                  xcart(:,list_symcoeff(2,icoeff2,1)),rprimd,&
2110 &                  cell(:,list_symcoeff(4,icoeff,1)),cell(:,1))>=cutoff.or.&
2111 &         distance_supercell(xcart(:,list_symcoeff(3,icoeff,1)),&
2112 &                  xcart(:,list_symcoeff(3,icoeff2,1)),rprimd,&
2113 &                  cell(:,list_symcoeff(4,icoeff,1)),&
2114 &                  cell(:,list_symcoeff(4,icoeff2,1)))>=cutoff)then
2115 !TEST_AM
2116            compatibleCoeffs(icoeff,icoeff2) = 0
2117            compatibleCoeffs(icoeff2,icoeff) = 0
2118 !TEST_AM
2119          end if
2120        end if
2121      end do
2122    end do
2123 
2124 
2125 !  Compute all the combination of coefficient up to the given order  (get the number)
2126    if(need_verbose)then
2127      write(message,'(1a)')' Compute the number of possible combinations'
2128      call wrtout(std_out,message,'COLL')
2129    end if
2130 
2131    ABI_ALLOCATE(list_coeff,(0))
2132    ABI_ALLOCATE(list_combination,(0,0))
2133    icoeff  = 1
2134    icoeff2 = 0
2135    call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2136 &                   list_coeff,list_combination,icoeff,max_power_strain,icoeff2,natom,ncoeff_sym,&
2137 &                   nstr_sym,icoeff,nrpt,nsym,1,power_disps(1),power_disps(2),symbols,nbody=option,&
2138 &                   compute=.false.,anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2139 &                   only_odd_power=need_only_odd_power,only_even_power=need_only_even_power)
2140    ncombination  = icoeff2
2141    ABI_DEALLOCATE(list_coeff)
2142    ABI_DEALLOCATE(list_combination)
2143 
2144 !  Compute all the combination of coefficient up to the given order
2145    if(need_verbose)then
2146      write(message,'(1a)')' Compute the combinations'
2147      call wrtout(std_out,message,'COLL')
2148    end if
2149 
2150    ABI_ALLOCATE(list_coeff,(0))
2151    ABI_ALLOCATE(list_combination_tmp,(power_disps(2),ncombination))
2152 
2153    icoeff  = 1
2154    icoeff2 = 0
2155    list_combination_tmp = 0
2156    call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,&
2157 &                   list_coeff,list_combination_tmp,icoeff,max_power_strain,icoeff2,natom,&
2158 &                   ncoeff_sym,nstr_sym,ncombination,nrpt,nsym,1,power_disps(1),power_disps(2),&
2159 &                   symbols,nbody=option,compute=.true.,&
2160 &                   anharmstr=need_anharmstr,spcoupling=need_spcoupling,&
2161 &                   only_odd_power=need_only_odd_power,only_even_power=need_only_even_power)
2162    ABI_DEALLOCATE(list_coeff)
2163    ABI_DEALLOCATE(compatibleCoeffs)
2164 
2165  else
2166    ABI_ALLOCATE(list_combination_tmp,(1,1))
2167  end if
2168 
2169  ABI_DEALLOCATE(xcart)
2170  ABI_DEALLOCATE(xred)
2171 
2172 !MPI
2173  if(need_verbose .and. nproc > 1)then
2174    write(message,'(1a)')' Distribute the combinations over the CPU'
2175    call wrtout(std_out,message,'COLL')
2176  end if
2177 
2178  call xmpi_bcast(ncombination, master, comm, ierr)
2179 
2180  ncoeff_alone = mod(ncombination,nproc)
2181  my_ncoeff = int(aint(real(ncombination,sp)/(nproc)))
2182 
2183  if(my_rank >= (nproc-ncoeff_alone)) then
2184    my_ncoeff = my_ncoeff  + 1
2185  end if
2186 
2187 !Set the buffsize for mpi scatterv
2188  ABI_ALLOCATE(buffsize,(nproc))
2189  ABI_ALLOCATE(buffdispl,(nproc))
2190  do ii = 1,nproc
2191    buffsize(ii) = int(aint(real(ncombination,sp)/(nproc))*power_disps(2))
2192    if(ii > (nproc-ncoeff_alone)) then
2193      buffsize(ii) = buffsize(ii) + power_disps(2)
2194    end if
2195  end do
2196 
2197  buffdispl(1) = 0
2198  do ii = 2,nproc
2199    buffdispl(ii) = buffdispl(ii-1) + buffsize(ii-1)
2200  end do
2201 
2202  ABI_ALLOCATE(list_combination,(power_disps(2),my_ncoeff))
2203  list_combination = 0
2204 
2205  my_size = my_ncoeff*power_disps(2)
2206  call xmpi_scatterv(list_combination_tmp,buffsize,buffdispl,list_combination,my_size,master,&
2207 &                   comm,ierr)
2208 
2209  ABI_DEALLOCATE(buffdispl)
2210  ABI_DEALLOCATE(buffsize)
2211  ABI_DEALLOCATE(list_combination_tmp)
2212 
2213  !Compute the coefficients from the list of combination
2214  if(need_verbose .and. nproc > 1)then
2215    write(message,'(1a)')' Compute the coefficients'
2216    call wrtout(std_out,message,'COLL')
2217  end if
2218  ABI_ALLOCATE(coeffs_tmp,(my_ncoeff))
2219  nterm      = nsym
2220  ndisp_max  = power_disps(2)
2221  ncoeff_max = my_ncoeff
2222  ABI_ALLOCATE(terms,(nterm))
2223  do ii=1,my_ncoeff
2224    call generateTermsFromList(cell,list_combination(:,ii),list_symcoeff,list_symstr,ncoeff_sym,&
2225 &                             ndisp_max,nrpt,nstr_sym,nsym,nterm,terms)
2226    call polynomial_coeff_init(one,nterm,coeffs_tmp(ii),terms(1:nterm),check=.true.)
2227 !  Free the terms array
2228    do iterm=1,nterm
2229      call polynomial_term_free(terms(iterm))
2230    end do
2231  end do
2232  ABI_DEALLOCATE(terms)
2233  ABI_DEALLOCATE(cell)
2234 
2235 
2236  ABI_DEALLOCATE(list_combination)
2237  ABI_DEALLOCATE(list_symcoeff)
2238  ABI_DEALLOCATE(list_symstr)
2239 
2240 !Final tranfert
2241 !1- Count the total number of coefficient
2242  ncoeff = 0
2243  do icoeff=1,ncoeff_max
2244    if (abs(coeffs_tmp(icoeff)%coefficient) >tol16) then
2245      ncoeff = ncoeff + 1
2246    end if
2247  end do
2248 
2249 !Get the total number of coefficients
2250 !ncoeff_max is the number of total coefficients before the symetries check
2251 !ncoeff_tot is the number of total coefficients after the symetries check
2252  ncoeff_tot = ncoeff!set the output
2253  call xmpi_sum(ncoeff_tot,comm,ierr)
2254  call xmpi_sum(ncoeff_max,comm,ierr)
2255 
2256 !Need to redistribute the coefficients over the CPU
2257 !Get the list with the number of coeff on each CPU
2258 !In order to be abble to compute the my_coeffindexes array which is for example:
2259 ! if CPU0 has 200  Coeff and CPU1 has 203 Coeff then
2260 ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403
2261  if(need_verbose .and. nproc > 1)then
2262    write(message,'(1a)')' Redistribute the coefficients over the CPU'
2263    call wrtout(std_out,message,'COLL')
2264  end if
2265 
2266  ABI_ALLOCATE(buffdispl,(nproc))
2267  buffdispl = 0
2268  buffdispl(my_rank+1) = my_ncoeff
2269  call xmpi_sum(buffdispl,comm,ierr)
2270  ABI_ALLOCATE(my_coeffindexes,(my_ncoeff))
2271  ABI_ALLOCATE(my_coefflist,(my_ncoeff))
2272  my_coeffindexes = 0
2273  my_coefflist = 0
2274  do icoeff=1,my_ncoeff
2275    my_coefflist(icoeff) = icoeff
2276    if(my_rank==0) then
2277      my_coeffindexes(icoeff) = icoeff
2278    else
2279      my_coeffindexes(icoeff) = sum(buffdispl(1:my_rank)) + icoeff
2280    end if
2281  end do
2282  ABI_DEALLOCATE(buffdispl)
2283 
2284 !Compute the new number of coefficient per CPU
2285  if(need_distributed) then
2286    ncoeff_alone = mod(ncoeff_tot,nproc)
2287    my_newncoeff = int(aint(real(ncoeff_tot,sp)/(nproc)))
2288    if(my_rank >= (nproc-ncoeff_alone)) then
2289      my_newncoeff = my_newncoeff  + 1
2290    end if
2291  else
2292    my_newncoeff = ncoeff_tot
2293  end if
2294 
2295  ncoeff = my_newncoeff ! Set the output
2296 
2297 !2:compute the number of coefficients and the list of the corresponding
2298 !  coefficients for each CPU.
2299  ABI_ALLOCATE(my_newcoeffindexes,(my_newncoeff))
2300  if(need_distributed) then
2301    do icoeff=1,my_newncoeff
2302      if(my_rank >= (nproc-ncoeff_alone))then
2303        my_newcoeffindexes(icoeff)=int(aint(real(ncoeff_tot,sp)/(nproc)))*(my_rank)+&
2304 &                              (my_rank - (nproc-ncoeff_alone)) + icoeff
2305      else
2306        my_newcoeffindexes(icoeff)=(my_newncoeff)*(my_rank)  + icoeff
2307      end if
2308    end do
2309  else
2310    do icoeff=1,my_newncoeff
2311      my_newcoeffindexes(icoeff) = icoeff
2312    end do
2313  end if
2314 
2315 !2- Transfer
2316  if(.not.need_distributed)then
2317    if(.not.allocated(coefficients))then
2318      ABI_DATATYPE_ALLOCATE(coefficients,(my_newncoeff))
2319    end if
2320  end if
2321  icoeff  = 0! icoeff is the current index in the total list of coefficients
2322  icoeff2 = 0! icoeff2 is the current index in the output coefficients array on each CPU
2323  icoeff3 = 0! icoeff3 is the current index in total new list of coefficients
2324  rank_to_send_save = 0
2325 
2326  do icoeff=1,ncoeff_max
2327 !  Need to send the rank with the chosen coefficient
2328    rank_to_send = 0
2329    my_icoeff = 0
2330    do ii=1,my_ncoeff
2331      if (my_coeffindexes(ii)==icoeff) then
2332        my_icoeff = ii
2333        if (abs(coeffs_tmp(my_icoeff)%coefficient) > tol16)then
2334          rank_to_send = my_rank
2335        else
2336          rank_to_send = -1
2337 !        Free the coefficient
2338          call polynomial_coeff_free(coeffs_tmp(ii))
2339        end if
2340        exit
2341      end if
2342    end do
2343    call xmpi_sum(rank_to_send, comm, ierr)
2344 !  This coefficient is not compute
2345    if (rank_to_send == -1) cycle
2346 
2347 !  increase icoeff3
2348    icoeff3 = icoeff3 + 1
2349 
2350 !  Find the receiver CPU
2351    rank_to_receive = 0
2352    do ii=1,my_newncoeff
2353      if (my_newcoeffindexes(ii)==icoeff3) then
2354        rank_to_receive = my_rank
2355      end if
2356    end do
2357    call xmpi_sum(rank_to_receive, comm, ierr)
2358 
2359    if(need_distributed.and.rank_to_send /= rank_to_send_save) then
2360      if(my_rank == rank_to_send_save)then
2361        ABI_DATATYPE_DEALLOCATE(coeffs_tmp)!Free memory if the current CPU has already distribute
2362                                           !all its own coefficients
2363      end if
2364      rank_to_send_save = rank_to_send
2365    end if
2366 
2367    if(need_distributed.and.my_rank == rank_to_receive)then
2368      if(.not.allocated(coefficients))then
2369        ABI_DATATYPE_ALLOCATE(coefficients,(my_newncoeff))
2370      end if
2371    end if
2372 
2373    if (need_distributed)then
2374      if(my_rank==rank_to_send)then
2375        if(any(my_newcoeffindexes(:)==icoeff3))then
2376          icoeff2 = icoeff2 + 1
2377 !        Get the name of this coefficient
2378          call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2379          call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2380 &                                  coeffs_tmp(my_icoeff)%terms,name=name,check=.false.)
2381        else
2382          call polynomial_coeff_MPIsend(coeffs_tmp(my_icoeff), icoeff, rank_to_receive, comm)
2383        end if
2384 !      Free the coefficient
2385        call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2386      else
2387        if(any(my_newcoeffindexes(:)==icoeff3))then
2388          icoeff2 = icoeff2 + 1
2389          call polynomial_coeff_MPIrecv(coefficients(icoeff2), icoeff, rank_to_send, comm)
2390          call polynomial_coeff_getName(name,coefficients(icoeff2),symbols,recompute=.TRUE.)
2391          call polynomial_coeff_SetName(name,coefficients(icoeff2))
2392        end if
2393      end if
2394    else
2395      icoeff2 = icoeff2 + 1
2396 !    Get the name of this coefficient
2397      if(my_rank==rank_to_send)then
2398        call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.)
2399        call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),&
2400 &                                 coeffs_tmp(my_icoeff)%terms,name=name,check=.false.)
2401 !      Free the coefficient
2402        call polynomial_coeff_free(coeffs_tmp(my_icoeff))
2403      end if
2404      call polynomial_coeff_broadcast(coefficients(icoeff2),rank_to_send, comm)
2405    end if
2406  end do
2407 
2408  ! if(iam_master)then
2409  !   filename = "terms_set.xml"
2410  !   call polynomial_coeff_writeXML(coefficients,my_newncoeff,filename=filename)
2411  ! end if
2412   
2413  if(need_verbose)then
2414    write(message,'(1x,I0,2a)') ncoeff_tot,' coefficients generated ',ch10
2415    call wrtout(ab_out,message,'COLL')
2416    call wrtout(std_out,message,'COLL')
2417  end if
2418 
2419 !Final deallocation
2420  ABI_DEALLOCATE(symbols)
2421  ABI_DEALLOCATE(my_coeffindexes)
2422  ABI_DEALLOCATE(my_newcoeffindexes)
2423  ABI_DEALLOCATE(my_coefflist)
2424  if(allocated(coeffs_tmp))then
2425    ABI_DATATYPE_DEALLOCATE(coeffs_tmp)
2426  end if
2427  
2428 end subroutine polynomial_coeff_getNorder

m_polynomial_coeff/polynomial_coeff_getOrder1 [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

3124 subroutine polynomial_coeff_getOrder1(cell,coeffs_out,list_symcoeff,&
3125 &                                     natom,ncoeff_out,ncoeff,nrpt,nsym,&
3126 &                                     symbols)
3127 
3128 
3129 !This section has been created automatically by the script Abilint (TD).
3130 !Do not modify the following lines by hand.
3131 #undef ABI_FUNC
3132 #define ABI_FUNC 'polynomial_coeff_getOrder1'
3133 !End of the abilint section
3134 
3135  implicit none
3136 
3137 !Arguments ------------------------------------
3138 !scalars
3139  integer,intent(in)  :: natom,ncoeff,nsym,nrpt
3140  integer,intent(out) :: ncoeff_out
3141 !arrays
3142  integer,intent(in) :: cell(3,nrpt)
3143  integer,intent(in) :: list_symcoeff(6,ncoeff,nsym)
3144  character(len=5),intent(in) :: symbols(natom)
3145  type(polynomial_coeff_type),allocatable,intent(inout) :: coeffs_out(:)
3146 !Local variables-------------------------------
3147 !scalar
3148  integer :: ia,ib,icoeff,icoeff_tmp,irpt,irpt_ref
3149  integer :: isym,iterm,mu,ncoeff_max,ndisp,nstrain,nterm_max
3150  real(dp):: coefficient,weight
3151 !arrays
3152  integer,allocatable :: atindx(:,:),cells(:,:,:),dir_int(:)
3153  integer,allocatable :: power_disps(:),power_strain(:),strain(:)
3154  character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/)
3155  character(len=200):: name
3156  character(len=500) :: message
3157  type(polynomial_term_type),dimension(:),allocatable :: terms
3158  type(polynomial_coeff_type),allocatable :: coeffs_tmp(:)
3159 !TEST_AM
3160 ! character(len=fnlen) :: filename
3161 !TEST_AM
3162 ! *************************************************************************
3163 
3164 !Initialisation of variables
3165  nterm_max  = nsym
3166  ncoeff_max = ncoeff
3167  ndisp = 1
3168  nstrain = 0
3169  ABI_ALLOCATE(coeffs_tmp,(ncoeff_max))
3170  ABI_ALLOCATE(terms,(nterm_max))
3171 
3172 
3173  icoeff_tmp = 0
3174  ABI_ALLOCATE(atindx,(2,ndisp))
3175  ABI_ALLOCATE(cells,(3,2,ndisp))
3176  ABI_ALLOCATE(dir_int,(ndisp))
3177  ABI_ALLOCATE(power_disps,(ndisp))
3178  ABI_ALLOCATE(power_strain,(nstrain))
3179  ABI_ALLOCATE(strain,(nstrain))
3180 
3181 !Found the ref cell
3182  irpt_ref = 1
3183  do irpt=1,nrpt
3184    if(all(cell(:,irpt)==0))then
3185      irpt_ref = irpt
3186      exit
3187    end if
3188  end do
3189 
3190  write(message,'(3a)') " Irreductible coefficient and associated atom 1, atom 2 and direction:",ch10,&
3191 &                     " for the 1st order"
3192  call wrtout(std_out,message,'COLL')
3193 
3194  do icoeff=1,ncoeff
3195 !  Reset counter
3196    iterm = 0
3197    coefficient = one
3198    do isym=1,nsym
3199      ndisp   = 1
3200      nstrain = 0
3201      mu   = list_symcoeff(1,icoeff,isym)
3202      ia   = list_symcoeff(2,icoeff,isym)
3203      ib   = list_symcoeff(3,icoeff,isym)
3204      irpt = list_symcoeff(4,icoeff,isym)
3205      weight = list_symcoeff(5,icoeff,isym)
3206 !    Fill First term arrays
3207      atindx(1,1) = ia; atindx(2,1) = ib;
3208      dir_int(1) = mu
3209      power_disps(1)   = 1
3210      cells(:,1,1) = (/0,0,0/)
3211      cells(:,2,1) = cell(:,irpt)
3212      iterm = iterm + 1
3213      call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(iterm),&
3214 &                              power_disps,power_strain,strain,weight,check=.true.)
3215    end do!end do sym
3216 
3217    if(iterm > 0)then
3218 !  increase coefficients and set it
3219      icoeff_tmp = icoeff_tmp + 1
3220      call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),terms(1:iterm),check=.true.)
3221    end if
3222 
3223 !  Deallocate the terms
3224    do iterm=1,nterm_max
3225      call polynomial_term_free(terms(iterm))
3226    end do
3227  end do!end do coeff_sym
3228 
3229  ABI_DEALLOCATE(terms)
3230  ABI_DEALLOCATE(atindx)
3231  ABI_DEALLOCATE(cells)
3232  ABI_DEALLOCATE(dir_int)
3233  ABI_DEALLOCATE(power_disps)
3234  ABI_DEALLOCATE(power_strain)
3235  ABI_DEALLOCATE(strain)
3236 
3237 !Count the number of terms
3238  ncoeff_out = 0
3239  do icoeff_tmp=1,ncoeff_max
3240    if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3241      ncoeff_out = ncoeff_out + 1
3242    end if
3243  end do
3244 
3245 !Transfer in the final array
3246  ABI_ALLOCATE(coeffs_out,(ncoeff_out))
3247  icoeff = 0
3248  do icoeff_tmp=1,ncoeff_max
3249    if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then
3250 !    Get the name of this coefficient
3251      call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.)
3252 !    Increase icoeff and fill the coeffs_out array
3253      icoeff = icoeff + 1
3254       call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,&
3255  &                               coeffs_out(icoeff),coeffs_tmp(icoeff_tmp)%terms,name=name)
3256 
3257      write(message,'(2a)')' ',trim(name)
3258      call wrtout(std_out,message,'COLL')
3259 
3260      do iterm = 1,coeffs_tmp(icoeff_tmp)%nterm
3261        write(message,'(a,I0,a,I0,2a)') '    Atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(1,1),&
3262 &                       ' and atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(2,1),&
3263 &                       ' in the direction ',mutodir(coeffs_tmp(icoeff_tmp)%terms(iterm)%direction(1))
3264        if(any(coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(:,2,1)/=0))then
3265          write(message,'(2a,I0,a,I0,a,I0,a)') trim(message),' in the cell ',&
3266 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(1,2,1),' ',&
3267 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(2,2,1),' ',&
3268 &                                       coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(3,2,1),'.'
3269        end if
3270        call wrtout(std_out,message,'COLL')
3271      end do
3272    end if
3273  end do
3274 
3275 !TEST_AM
3276 ! filename = "terms_1st_order.xml"
3277 ! call polynomial_coeff_writeXML(coeffs_out,ncoeff_out,filename=filename)
3278 !TEST_AM
3279 
3280  write(message,'(a,1x,I0,a)') ch10,&
3281 &       ncoeff_out,' fitted coefficients for the 1st order '
3282  call wrtout(ab_out,message,'COLL')
3283  call wrtout(std_out,message,'COLL')
3284 
3285 !Deallocation
3286  do icoeff=1,ncoeff_max
3287    call polynomial_coeff_free(coeffs_tmp(icoeff))
3288  end do
3289  ABI_DEALLOCATE(coeffs_tmp)
3290 
3291 end subroutine polynomial_coeff_getOrder1

m_polynomial_coeff/polynomial_coeff_init [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_anharmonics_terms,m_effective_potential_file,m_fit_polynomial_coeff
      m_polynomial_coeff,mover_effpot

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

132 subroutine polynomial_coeff_init(coefficient,nterm,polynomial_coeff,terms,name,check)
133 
134 
135 !This section has been created automatically by the script Abilint (TD).
136 !Do not modify the following lines by hand.
137 #undef ABI_FUNC
138 #define ABI_FUNC 'polynomial_coeff_init'
139 !End of the abilint section
140 
141  implicit none
142 
143 !Arguments ------------------------------------
144 !scalars
145  integer, intent(in) :: nterm
146  real(dp),intent(in) :: coefficient
147  logical,optional,intent(in) :: check
148 !arrays
149  character(len=200),optional,intent(in) :: name
150  type(polynomial_term_type),intent(in) :: terms(nterm)
151  type(polynomial_coeff_type), intent(out) :: polynomial_coeff
152 !Local variables-------------------------------
153 !scalar
154  integer :: iterm1,iterm2
155  integer :: ii,nterm_tmp
156  real(dp):: coefficient_tmp
157  logical :: check_in = .false.
158 !arrays
159  real(dp) :: weights(nterm)
160  character(len=200) :: name_tmp
161 ! *************************************************************************
162 
163 !First free before initilisation
164  call polynomial_coeff_free(polynomial_coeff)
165 
166  if(present(check)) check_in = check
167  if(check_in)then
168 !  Check if the list of term is available or contains identical terms
169 !  in this case, remove all the not needed terms
170    nterm_tmp = 0
171    weights(:) = one
172    do iterm1=1,nterm
173      if(abs(weights(iterm1)) < tol16)cycle
174      weights(iterm1) = terms(iterm1)%weight
175      do iterm2=iterm1+1,nterm
176        if(abs(weights(iterm2)) < tol16)cycle
177 !      if the terms are identical we check the weight
178        if(terms(iterm1)==terms(iterm2))then
179          weights(iterm1) = weights(iterm1) + terms(iterm2)%weight
180          weights(iterm2) = 0
181        end if
182      end do
183      if(abs(weights(iterm1)) > tol16) then
184        weights(iterm1)= anint(weights(iterm1)/weights(iterm1))
185      end if
186    end do
187 
188 !  Count the number of terms
189    nterm_tmp = 0
190    do iterm1=1,nterm
191      if(abs(weights(iterm1)) > tol16)then
192        nterm_tmp = nterm_tmp + 1
193      end if
194    end do
195 
196    if (nterm_tmp ==0)then
197      coefficient_tmp = 0.0
198    else
199      coefficient_tmp = coefficient
200    end if
201 
202  else
203    nterm_tmp = nterm
204    coefficient_tmp = coefficient
205    weights(:) = terms(:)%weight
206  end if!end Check
207 
208  if(present(name))then
209    name_tmp = name
210  else
211    name_tmp = ""
212  end if
213 
214 !Initilisation
215  polynomial_coeff%name = name_tmp
216  polynomial_coeff%nterm = nterm_tmp
217  polynomial_coeff%coefficient = coefficient_tmp
218  ABI_DATATYPE_ALLOCATE(polynomial_coeff%terms,(polynomial_coeff%nterm))
219  iterm1 = 0
220  do ii = 1,nterm
221    if(abs(weights(ii)) > tol16)then
222      iterm1 = iterm1 + 1
223      call polynomial_term_init(terms(ii)%atindx,terms(ii)%cell,terms(ii)%direction,terms(ii)%ndisp,&
224 &                              terms(ii)%nstrain,polynomial_coeff%terms(iterm1),terms(ii)%power_disp,&
225 &                              terms(ii)%power_strain,terms(ii)%strain,terms(ii)%weight)
226    end if
227  end do
228 
229 end subroutine polynomial_coeff_init

m_polynomial_coeff/polynomial_coeff_MPIrecv [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

741 subroutine polynomial_coeff_MPIrecv(coefficients, tag, source, comm)
742 
743 
744 !This section has been created automatically by the script Abilint (TD).
745 !Do not modify the following lines by hand.
746 #undef ABI_FUNC
747 #define ABI_FUNC 'polynomial_coeff_MPIrecv'
748 !End of the abilint section
749 
750  implicit none
751 
752 !Arguments ------------------------------------
753 !array
754  type(polynomial_coeff_type),intent(inout) :: coefficients
755  integer, intent(in) :: source,comm,tag
756 
757 !Local variables-------------------------------
758 !scalars
759  integer :: ierr,ii
760 !arrays
761 
762 ! *************************************************************************
763 
764  if (xmpi_comm_size(comm) == 1) return
765 
766 
767 ! Free the output
768   call polynomial_coeff_free(coefficients)
769 
770  ! Transmit variables
771   call xmpi_recv(coefficients%name, source, 9*tag+0, comm, ierr)
772   call xmpi_recv(coefficients%nterm, source, 9*tag+1, comm, ierr)
773   call xmpi_recv(coefficients%coefficient, source, 9*tag+2, comm, ierr)
774 
775  !Allocate arrays on the other nodes.
776   ABI_DATATYPE_ALLOCATE(coefficients%terms,(coefficients%nterm))
777   do ii=1,coefficients%nterm
778     call polynomial_term_free(coefficients%terms(ii))
779   end do
780 
781 ! Set the number of term on each node (needed for allocations of array)
782   do ii = 1,coefficients%nterm
783     call xmpi_recv(coefficients%terms(ii)%ndisp, source, 9*tag+3, comm, ierr)
784     call xmpi_recv(coefficients%terms(ii)%nstrain, source, 9*tag+4, comm, ierr)
785   end do
786 
787 ! Allocate arrays on the other nodes
788   do ii = 1,coefficients%nterm
789     ABI_ALLOCATE(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp))
790     coefficients%terms(ii)%atindx = 0
791     ABI_ALLOCATE(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp))
792     ABI_ALLOCATE(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp))
793     ABI_ALLOCATE(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp))
794     ABI_ALLOCATE(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain))
795     ABI_ALLOCATE(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain))
796   end do
797 
798 ! Transfert value
799   do ii = 1,coefficients%nterm
800     call xmpi_recv(coefficients%terms(ii)%weight, source, 9*tag+5, comm, ierr)
801     call xmpi_recv(coefficients%terms(ii)%atindx, source, 9*tag+6, comm, ierr)
802     call xmpi_recv(coefficients%terms(ii)%direction, source, 9*tag+7, comm, ierr)
803     call xmpi_recv(coefficients%terms(ii)%cell, source, 9*tag+8, comm, ierr)
804     call xmpi_recv(coefficients%terms(ii)%power_disp, source, 9*tag+9, comm, ierr)
805     call xmpi_recv(coefficients%terms(ii)%power_strain, source, 9*tag+10, comm, ierr)
806     call xmpi_recv(coefficients%terms(ii)%strain, source, 9*tag+11, comm, ierr)
807   end do
808 
809 end subroutine polynomial_coeff_MPIrecv

m_polynomial_coeff/polynomial_coeff_MPIsend [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

665 subroutine polynomial_coeff_MPIsend(coefficients, tag, dest, comm)
666 
667 
668 !This section has been created automatically by the script Abilint (TD).
669 !Do not modify the following lines by hand.
670 #undef ABI_FUNC
671 #define ABI_FUNC 'polynomial_coeff_MPIsend'
672 !End of the abilint section
673 
674  implicit none
675 
676 !Arguments ------------------------------------
677 !array
678  type(polynomial_coeff_type),intent(inout) :: coefficients
679  integer, intent(in) :: dest,comm,tag
680 
681 !Local variables-------------------------------
682 !scalars
683  integer :: ierr,ii
684  integer :: my_rank
685 !arrays
686 
687 ! *************************************************************************
688 
689  if (xmpi_comm_size(comm) == 1) return
690 
691   my_rank = xmpi_comm_rank(comm)
692 ! Transmit variables
693   call xmpi_send(coefficients%name, dest, 9*tag+0, comm, ierr)
694   call xmpi_send(coefficients%nterm, dest, 9*tag+1, comm, ierr)
695   call xmpi_send(coefficients%coefficient, dest, 9*tag+2, comm, ierr)
696 
697 ! Set the number of term on each node (needed for allocations of array)
698   do ii = 1,coefficients%nterm
699     call xmpi_send(coefficients%terms(ii)%ndisp, dest, 9*tag+3, comm, ierr)
700     call xmpi_send(coefficients%terms(ii)%nstrain, dest, 9*tag+4, comm, ierr)
701   end do
702 
703 ! Transfert value
704   do ii = 1,coefficients%nterm
705       call xmpi_send(coefficients%terms(ii)%weight, dest, 9*tag+5, comm, ierr)
706       call xmpi_send(coefficients%terms(ii)%atindx, dest, 9*tag+6, comm, ierr)
707       call xmpi_send(coefficients%terms(ii)%direction, dest, 9*tag+7, comm, ierr)
708       call xmpi_send(coefficients%terms(ii)%cell, dest, 9*tag+8, comm, ierr)
709       call xmpi_send(coefficients%terms(ii)%power_disp, dest, 9*tag+9, comm, ierr)
710       call xmpi_send(coefficients%terms(ii)%power_strain, dest, 9*tag+10, comm, ierr)
711       call xmpi_send(coefficients%terms(ii)%strain, dest, 9*tag+11, comm, ierr)
712   end do
713 
714 end subroutine polynomial_coeff_MPIsend

m_polynomial_coeff/polynomial_coeff_setCoefficient [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_effective_potential_file,mover_effpot

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

314 subroutine polynomial_coeff_setCoefficient(coefficient,polynomial_coeff)
315 
316 
317 !This section has been created automatically by the script Abilint (TD).
318 !Do not modify the following lines by hand.
319 #undef ABI_FUNC
320 #define ABI_FUNC 'polynomial_coeff_setCoefficient'
321 !End of the abilint section
322 
323  implicit none
324 
325 !Arguments ------------------------------------
326 !scalars
327  real(dp),intent(in) :: coefficient
328 !arrays
329  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
330 !Local variables-------------------------------
331 !scalar
332 !arrays
333 ! *************************************************************************
334 
335  polynomial_coeff%coefficient = coefficient
336 
337 end subroutine polynomial_coeff_setCoefficient

m_polynomial_coeff/polynomial_coeff_setName [ Functions ]

[ Top ] [ m_polynomial_coeff ] [ 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

PARENTS

      m_effective_potential_file,m_fit_polynomial_coeff,m_polynomial_coeff

CHILDREN

      polynomial_coeff_free,polynomial_coeff_getname,polynomial_coeff_init
      polynomial_term_free,polynomial_term_init,wrtout

SOURCE

363 subroutine polynomial_coeff_setName(name,polynomial_coeff)
364 
365 
366 !This section has been created automatically by the script Abilint (TD).
367 !Do not modify the following lines by hand.
368 #undef ABI_FUNC
369 #define ABI_FUNC 'polynomial_coeff_setName'
370 !End of the abilint section
371 
372  implicit none
373 
374 !Arguments ------------------------------------
375 !scalars
376 !arrays
377  character(len=200),intent(in) :: name
378  type(polynomial_coeff_type), intent(inout) :: polynomial_coeff
379 !Local variables-------------------------------
380 !scalar
381 !arrays
382 ! *************************************************************************
383 
384  polynomial_coeff%name = name
385 
386 end subroutine polynomial_coeff_setName

m_polynomial_coeff/polynomial_coeff_type [ Types ]

[ Top ] [ m_polynomial_coeff ] [ 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

77  type, public :: polynomial_coeff_type
78 
79    character(len=200) :: name = ""
80 !     Name of the polynomial_coeff (Sr_y-O1_y)^3) for example
81 
82    integer :: nterm = 0
83 !     Number of terms (short range interaction) for this polynomial_coeff
84 
85    real(dp) :: coefficient = zero
86 !     coefficient = value of the coefficient of this term
87 !     \frac{\partial E^{k}}{\partial \tau^{k}}
88 
89    type(polynomial_term_type),dimension(:),allocatable :: terms
90 !     polynomial_term(nterm)<type(polynomial_term)>
91 !     contains all the displacements for this coefficient
92 
93  end type polynomial_coeff_type