TABLE OF CONTENTS


ABINIT/m_effective_potential_file [ Modules ]

[ Top ] [ Modules ]

NAME

 m_effective_potential_file

FUNCTION

 This  module contains all routine to read the effective potential from files
 Can also read coefficients from XML
 (XML or DDB)

COPYRIGHT

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

PARENTS

SOURCE

 22 #if defined HAVE_CONFIG_H
 23 #include "config.h"
 24 #endif
 25 
 26 #include "abi_common.h"
 27 
 28 module m_effective_potential_file
 29 
 30  use defs_basis
 31  use m_errors
 32  use m_abicore
 33  use m_xmpi
 34  use m_harmonics_terms
 35  use m_anharmonics_terms
 36  use m_effective_potential
 37  use m_ifc
 38 #if defined HAVE_NETCDF
 39  use netcdf
 40 #endif
 41 
 42  use m_io_tools,   only : open_file
 43  use m_geometry,   only : xcart2xred, metric
 44  use m_symfind,    only : symfind, symlatt
 45  use m_crystal,    only : crystal_t, crystal_init, crystal_free
 46  use m_dynmat,     only : dfpt_prtph
 47  use m_abihist,    only : abihist,abihist_init,abihist_free,abihist_copy,read_md_hist
 48  use m_ddb_internalstr, only : ddb_internalstr
 49 
 50  implicit none
 51 
 52  public :: effective_potential_file_getDimCoeff
 53  public :: effective_potential_file_getDimMD
 54  public :: effective_potential_file_getDimSystem
 55  public :: effective_potential_file_getDimStrainCoupling
 56  public :: effective_potential_file_getType
 57  public :: effective_potential_file_mapHistToRef
 58  public :: effective_potential_file_read
 59  public :: effective_potential_file_readDisplacement
 60  public :: effective_potential_file_readMDfile
 61  private :: coeffs_xml2effpot
 62  private :: system_getDimFromXML
 63  private :: system_xml2effpot
 64  private :: system_ddb2effpot
 65 
 66 #ifndef HAVE_LIBXML
 67  private :: rdfromline
 68  private :: rmtabfromline
 69  private :: rdfromline_value
 70  private :: elementfromline
 71 #endif
 72 
 73 #if defined HAVE_LIBXML
 74  public :: effpot_xml_checkXML
 75  public :: effpot_xml_getDimCoeff
 76  public :: effpot_xml_readSystem
 77  public :: effpot_xml_getValue
 78  public :: effpot_xml_getAttribute
 79  public :: effpot_xml_getDimSystem
 80 
 81  interface
 82    subroutine effpot_xml_readSystem(filename,natom,&
 83 &     ntypat,nrpt,nqpt,amu,atmfrc,cell,dynmat,elastic_constants,&
 84 &     energy,epsilon_inf,ewald_atmfrc,&
 85 &     phfrq,rprimd,qph1l,short_atmfrc,typat,xcart,zeff)&
 86 &                          bind(C,name="effpot_xml_readSystem")
 87      use iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
 88      integer(C_INT) :: natom,ntypat,nrpt,nqpt
 89      integer(C_INT) :: typat(natom)
 90      integer(C_INT) :: cell(3,nrpt)
 91      real(C_DOUBLE) :: energy
 92      real(C_DOUBLE) :: dynmat(2,3,natom,3,natom,nqpt)
 93      real(C_DOUBLE) :: phfrq(3*natom,nqpt),qph1l(3,nqpt)
 94      real(C_DOUBLE) :: atmfrc(3,natom,3,natom,nrpt)
 95      real(C_DOUBLE) :: short_atmfrc(3,natom,3,natom,nrpt)
 96      real(C_DOUBLE) :: ewald_atmfrc(3,natom,3,natom,nrpt)
 97      real(C_DOUBLE) :: amu(ntypat),rprimd(3,3),epsilon_inf(3,3)
 98      real(C_DOUBLE) :: zeff(3,3,natom)
 99      real(C_DOUBLE) :: elastic_constants(6,6),xcart(3,natom)
100      character(kind=C_CHAR) :: filename(*)
101    end subroutine effpot_xml_readSystem
102  end interface
103 
104  interface
105    subroutine effpot_xml_readStrainCoupling(filename,natom,&
106 &     nrpt,voigt,elastic3rd,elastic_displacement,&
107 &     strain_coupling,phonon_strain_atmfrc,phonon_straincell)&
108 &                          bind(C,name="effpot_xml_readStrainCoupling")
109      use iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
110      integer(C_INT) :: natom
111      integer(C_INT) :: nrpt,voigt
112      integer(c_INT) :: phonon_straincell(3,nrpt)
113      real(C_DOUBLE) :: elastic3rd(6,6),elastic_displacement(6,3,natom)
114      real(C_DOUBLE) :: strain_coupling(3,natom)
115      real(C_DOUBLE) :: phonon_strain_atmfrc(3,natom,3,natom,nrpt)
116      character(kind=C_CHAR) :: filename(*)
117    end subroutine effpot_xml_readStrainCoupling
118  end interface
119 
120  interface
121    subroutine effpot_xml_readCoeff(filename,ncoeff,ndisp,nterm,&
122 &                                 coefficient,atindx,cell,direction,power_disp,&
123 &                                 power_strain,strain,weight)&
124 &                          bind(C,name="effpot_xml_readCoeff")
125      use iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT
126      character(kind=C_CHAR) :: filename(*)
127      integer(C_INT) :: ncoeff,ndisp,nterm
128      integer(C_INT) :: atindx(ncoeff,nterm,2,ndisp)
129      integer(C_INT) :: cell(ncoeff,nterm,3,2,ndisp)
130      integer(C_INT) :: direction(ncoeff,nterm,ndisp)
131      integer(C_INT) :: strain(ncoeff,nterm,ndisp)
132      integer(C_INT) :: power_disp(ncoeff,nterm,ndisp)
133      integer(C_INT) :: power_strain(ncoeff,nterm,ndisp)
134      real(C_DOUBLE) :: coefficient(ncoeff)
135      real(C_DOUBLE) :: weight(ncoeff,nterm)
136    end subroutine effpot_xml_readCoeff
137  end interface
138 
139  interface
140    subroutine effpot_xml_getDimSystem(filename,natom,ntypat,nqpt,nrpt1,nrpt2)&
141 &                          bind(C,name="effpot_xml_getDimSystem")
142      use iso_c_binding, only : C_CHAR,C_INT
143      integer(C_INT) :: natom,ntypat,nqpt,nrpt1,nrpt2
144      character(kind=C_CHAR) :: filename(*)
145    end subroutine effpot_xml_getDimSystem
146  end interface
147 
148  interface
149    subroutine effpot_xml_getDimStrainCoupling(filename,nrpt,voigt)&
150 &                          bind(C,name="effpot_xml_getDimStrainCoupling")
151      use iso_c_binding, only : C_CHAR,C_INT
152      integer(C_INT) :: voigt
153      integer(C_INT) :: nrpt
154      character(kind=C_CHAR) :: filename(*)
155    end subroutine effpot_xml_getDimStrainCoupling
156  end interface
157 
158  interface
159    subroutine effpot_xml_getDimCoeff(filename,ncoeff,nterm_max,ndisp_max)&
160 &                          bind(C,name="effpot_xml_getDimCoeff")
161      use iso_c_binding, only : C_CHAR,C_DOUBLE,C_INT,C_PTR
162      character(kind=C_CHAR) :: filename(*)
163 !     character(kind=C_CHAR) :: name(*)
164      type(C_PTR) :: name
165      integer(C_INT) :: ncoeff,ndisp_max,nterm_max
166    end subroutine effpot_xml_getDimCoeff
167  end interface
168 
169 
170  interface
171    subroutine effpot_xml_checkXML(filename,name_root) &
172 &                          bind(C,name="effpot_xml_checkXML")
173      use iso_c_binding, only : C_CHAR
174      character(kind=C_CHAR) :: filename(*),name_root(*)
175    end subroutine effpot_xml_checkXML
176  end interface
177 
178  interface
179    subroutine effpot_xml_getValue(filename,name_value,value_result) &
180  &                          bind(C,name="effpot_xml_getValue")
181       use iso_c_binding, only : C_CHAR
182       implicit none
183       character(kind=C_CHAR) :: filename(*),name_value(*)
184       character(kind=C_CHAR) :: value_result
185     end subroutine effpot_xml_getValue
186   end interface
187 
188  interface
189    subroutine effpot_xml_getAttribute(filename,name_value,name_attribute) &
190 &                          bind(C,name="effpot_xml_getAttribute")
191      use iso_c_binding, only : C_CHAR
192      character(kind=C_CHAR) :: filename(*),name_value(*),name_attribute(*)
193    end subroutine effpot_xml_getAttribute
194  end interface
195 
196  interface
197    subroutine effpot_xml_getNumberKey(filename,name_value,number) &
198 &                          bind(C,name="effpot_xml_getNumberKey")
199      use iso_c_binding, only : C_CHAR,C_INT
200      character(kind=C_CHAR) :: filename(*),name_value(*)
201      integer(C_INT) :: number
202    end subroutine effpot_xml_getNumberKey
203  end interface
204 
205 #endif
206 
207  integer,parameter :: XML_RECL = 50000

m_effective_potential_file/coeffs_xml2effpot [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 coeffs_xml2effpot

FUNCTION

 Open xml file of effective potentiel, then reads the variables
 and store them in effective potentential type

INPUTS

 filename = path of input or output file
 comm=MPI communicator

OUTPUT

 eff_pot<type(effective_potential_type)> = effective_potential datatype

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

2977 subroutine coeffs_xml2effpot(eff_pot,filename,comm)
2978 
2979  use m_atomdata
2980  use m_effective_potential, only : effective_potential_type
2981  use m_polynomial_coeff
2982  use m_polynomial_term
2983  use m_crystal, only : symbols_crystal
2984 #if defined HAVE_LIBXML
2985  use iso_c_binding, only : C_CHAR,C_PTR,c_f_pointer
2986 #endif
2987 
2988 !This section has been created automatically by the script Abilint (TD).
2989 !Do not modify the following lines by hand.
2990 #undef ABI_FUNC
2991 #define ABI_FUNC 'coeffs_xml2effpot'
2992 !End of the abilint section
2993 
2994  implicit none
2995 
2996  !Arguments ------------------------------------
2997  !scalars
2998  character(len=*),intent(in) :: filename
2999  integer, intent(in) :: comm
3000  !arrays
3001  type(effective_potential_type), intent(inout) :: eff_pot
3002 
3003  !Local variables-------------------------------
3004  !scalar
3005  integer :: ii,jj,my_rank,ndisp,ncoeff,nterm_max,nstrain,ndisp_max,nproc,nterm
3006 ! character(len=200),allocatable :: name(:)
3007  character(len=200) :: name
3008 #ifdef HAVE_LIBXML
3009  integer :: icoeff,iterm
3010 #endif
3011 
3012 #ifndef HAVE_LIBXML
3013  integer :: funit = 1,ios = 0
3014  integer :: icoeff,idisp,istrain,iterm,mu
3015  logical :: found,found2,displacement
3016  character (len=XML_RECL) :: line,readline
3017  character (len=XML_RECL) :: strg,strg1
3018 #endif
3019  character(len=500) :: message
3020  character(len=5),allocatable :: symbols(:)
3021  integer,parameter :: master=0
3022  logical :: iam_master
3023  logical :: debug =.FALSE.
3024  !arrays
3025  real(dp),allocatable :: coefficient(:),weight(:,:)
3026  integer,allocatable :: atindx(:,:,:,:), cell(:,:,:,:,:),direction(:,:,:),power_disp(:,:,:)
3027  integer,allocatable :: strain(:,:,:),power_strain(:,:,:)
3028  type(polynomial_coeff_type),dimension(:),allocatable :: coeffs
3029  type(polynomial_term_type),dimension(:,:),allocatable :: terms
3030 
3031 ! *************************************************************************
3032 
3033 
3034  !Open the atomicdata XML file for reading
3035  write(message,'(a,a)')'-Opening the file ',filename
3036 
3037  call wrtout(ab_out,message,'COLL')
3038  call wrtout(std_out,message,'COLL')
3039 
3040  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3041  iam_master = (my_rank == master)
3042 
3043 !Get Dimention of system and allocation/initialisation of array
3044  ncoeff  = 0
3045  nterm   = 0
3046  ndisp   = 0
3047  nstrain = 0
3048  call effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max)
3049 
3050 !  Do some checks
3051  if (nterm_max<=0) then
3052    write(message, '(a,a,a)' )&
3053 &     ' Unable to read the number of terms in ',trim(filename),ch10
3054    MSG_ERROR(message)
3055  end if
3056 
3057   if (ndisp_max<=0) then
3058     write(message, '(a,a,a)' )&
3059 &    ' Unable to read the number of displacement in ',trim(filename),ch10
3060     MSG_ERROR(message)
3061   end if
3062 
3063 !Allocation ov the polynomial coeff type
3064  ABI_DATATYPE_ALLOCATE(coeffs,(ncoeff))
3065 
3066  if(iam_master)then
3067 
3068 #if defined HAVE_LIBXML
3069    write(message,'(3a)')'-Reading the file ',trim(filename),&
3070 &   ' with LibXML library'
3071 #else
3072    write(message,'(3a)')'-Reading the file ',trim(filename),&
3073 &   ' with Fortran'
3074 #endif
3075    call wrtout(ab_out,message,'COLL')
3076    call wrtout(std_out,message,'COLL')
3077 
3078 
3079    ABI_ALLOCATE(symbols,(eff_pot%crystal%natom))
3080 !  Get the symbols arrays
3081    call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,&
3082 &                       eff_pot%crystal%npsp,symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
3083 
3084 
3085  !Read with libxml librarie
3086 #if defined HAVE_LIBXML
3087 
3088    ABI_DATATYPE_ALLOCATE(terms,(ncoeff,nterm_max))
3089    ABI_ALLOCATE(atindx,(ncoeff,nterm_max,2,ndisp_max))
3090    ABI_ALLOCATE(coefficient,(ncoeff))
3091    ABI_ALLOCATE(cell,(ncoeff,nterm_max,3,2,ndisp_max))
3092    ABI_ALLOCATE(direction,(ncoeff,nterm_max,ndisp_max))
3093    ABI_ALLOCATE(strain,(ncoeff,nterm_max,ndisp_max))
3094    ABI_ALLOCATE(power_disp,(ncoeff,nterm_max,ndisp_max))
3095    ABI_ALLOCATE(power_strain,(ncoeff,nterm_max,ndisp_max))
3096    ABI_ALLOCATE(weight,(ncoeff,nterm_max))
3097 
3098 !  Read the values of this term with libxml
3099    call effpot_xml_readCoeff(char_f2c(trim(filename)),ncoeff,ndisp_max,nterm_max,&
3100 &                            coefficient,atindx,cell,direction,power_disp,power_strain,&
3101 &                            strain,weight)
3102 !  In the XML the atom index begin to zero
3103 !  Need to shift for fortran array
3104    atindx(:,:,:,:) = atindx(:,:,:,:) + 1
3105 
3106    do icoeff=1,ncoeff
3107      do iterm=1,nterm_max
3108 !      Initialisation of the polynomial_term structure with the values from the
3109        call polynomial_term_init(atindx(icoeff,iterm,:,:),cell(icoeff,iterm,:,:,:),&
3110 &                                direction(icoeff,iterm,:),ndisp_max,ndisp_max,terms(icoeff,iterm),&
3111 &                                power_disp(icoeff,iterm,:),power_strain(icoeff,iterm,:),&
3112 &                                strain(icoeff,iterm,:),weight(icoeff,iterm),check=.true.)
3113      end do
3114 !    Initialisation of the polynomial_coefficent structure with the values
3115      call polynomial_coeff_init(coefficient(icoeff),nterm_max,coeffs(icoeff),&
3116 &                               terms(icoeff,:),check=.true.)
3117 
3118 !    Get the name of this coefficient  and set it
3119 !    Try to find the index of the term corresponding to the interation in the
3120 !    reference cell (000) in order to compute the name correctly...
3121 !    If this coeff is not in the ref cell, take by default the first term
3122      if(coeffs(icoeff)%nterm > 0)then
3123        call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.)
3124        call polynomial_coeff_setName(name,coeffs(icoeff))
3125      end if
3126 
3127 !    Free them all
3128      do iterm=1,nterm_max
3129        call polynomial_term_free(terms(icoeff,iterm))
3130      end do
3131    end do
3132 
3133 #else
3134    ABI_DATATYPE_ALLOCATE(terms,(1,nterm_max))
3135    ABI_ALLOCATE(atindx,(1,1,2,ndisp_max))
3136    ABI_ALLOCATE(coefficient,(1))
3137    ABI_ALLOCATE(cell,(1,1,3,2,ndisp_max))
3138    ABI_ALLOCATE(direction,(1,1,ndisp_max))
3139    ABI_ALLOCATE(strain,(1,1,ndisp_max))
3140    ABI_ALLOCATE(power_disp,(1,1,ndisp_max))
3141    ABI_ALLOCATE(power_strain,(1,1,ndisp_max))
3142    ABI_ALLOCATE(weight,(1,1))
3143 !  Loop over the file
3144 !  Read the values of all the terms with fortran
3145    if (open_file(filename,message,unit=funit,form="formatted",&
3146 &              status="old",action="read") /= 0) then
3147      MSG_ERROR(message)
3148    end if
3149 
3150 !    Start a reading loop in fortran
3151      rewind(unit=funit)
3152      ios  = 0
3153      found=.false.
3154 
3155 !    Initialisation of counter
3156      icoeff  = 0
3157 
3158 !    Parser
3159      do while (ios==0)
3160        read(funit,'(a)',iostat=ios) readline
3161        if (ios == 0) then
3162          call rmtabfromline(readline)
3163          line=adjustl(readline)
3164          if ((line(1:12)==char(60)//'coefficient')) then
3165 !          Read headers of coefficient
3166            call rdfromline('text',line,strg)
3167            if (strg/="") then
3168              name=trim(strg)
3169            end if
3170            call rdfromline('value',line,strg)
3171            if (strg/="") then
3172              strg1=trim(strg)
3173              read(strg1,*) coefficient(1)
3174            else
3175              coefficient(1) = zero
3176            end if
3177 !          End read headers of coefficient
3178 !          Reset counter
3179            found  = .false.
3180            atindx = 0;  cell   = 0 ;  direction = 0
3181            strain = 0; power_strain = 0;  power_disp  = 0
3182            iterm   = 0
3183            idisp   = 0
3184            istrain = 0
3185            nterm   = 0
3186            do while (.not.found)
3187              read(funit,'(a)',iostat=ios) readline
3188              call rmtabfromline(readline)
3189              line=adjustl(readline)
3190              if ((line(1:13)==char(60)//'/coefficient')) then
3191                found= .true.
3192                cycle
3193              end if
3194              if ((line(1:5)==char(60)//'term')) then
3195                nterm = nterm + 1
3196                ndisp = 0
3197                nstrain = 0
3198                idisp = 0
3199                istrain = 0
3200                displacement = .true.
3201                call rdfromline('weight',line,strg)
3202                if (strg/="") then
3203                  strg1=trim(strg)
3204                  read(strg1,*) weight
3205                end if
3206                do while(displacement)
3207                  read(funit,'(a)',iostat=ios) readline
3208                  call rmtabfromline(readline)
3209                  line=adjustl(readline)
3210                  if ((line(1:6)==char(60)//'/term')) then
3211                    displacement = .false.
3212                  end if
3213                  if ((line(1:7)==char(60)//'strain')) then
3214                    nstrain = nstrain + 1
3215                    istrain = istrain + 1
3216                    call rdfromline('power',line,strg)
3217                    if (strg/="") then
3218                      strg1=trim(strg)
3219                      read(strg1,*) power_strain(1,1,istrain)
3220                    end if
3221                    call rdfromline('voigt',line,strg)
3222                    if (strg/="") then
3223                      strg1=trim(strg)
3224                      read(strg1,*) strain(1,1,istrain)
3225                    end if
3226                  end if
3227                  if ((line(1:18)==char(60)//'displacement_diff')) then
3228                    ndisp = ndisp + 1
3229                    idisp = idisp + 1
3230                    found2=.true.
3231                    call rdfromline('atom_a',line,strg)
3232                    if (strg/="") then
3233                      strg1=trim(strg)
3234                      read(strg1,*) atindx(1,1,1,idisp)
3235                    end if
3236                    call rdfromline('atom_b',line,strg)
3237                    if (strg/="") then
3238                      strg1=trim(strg)
3239                      read(strg1,*) atindx(1,1,2,idisp)
3240                    end if
3241                    call rdfromline('direction',line,strg)
3242                    if (strg/="") then
3243                      strg1=trim(strg)
3244                      if (trim(strg1).eq."x") direction(1,1,idisp) = 1
3245                      if (trim(strg1).eq."y") direction(1,1,idisp) = 2
3246                      if (trim(strg1).eq."z") direction(1,1,idisp) = 3
3247                    end if
3248                    call rdfromline('power',line,strg)
3249                    if (strg/="") then
3250                      strg1=trim(strg)
3251                      read(strg1,*) power_disp(1,1,idisp)
3252                    end if
3253                    do while(found2)
3254                      read(funit,'(a)',iostat=ios) readline
3255                      call rmtabfromline(readline)
3256                      line=adjustl(readline)
3257                      if ((line(1:7)==char(60)//'cell_a')) then
3258                        call rdfromline_value('cell_a',line,strg)
3259                        if (strg/="") then
3260                          strg1=trim(strg)
3261                          read(strg1,*) (cell(1,1,mu,1,idisp),mu=1,3)
3262                        else
3263                          read(funit,'(a)',iostat=ios) readline
3264                          call rmtabfromline(readline)
3265                          line=adjustl(readline)
3266                          call rdfromline_value('cell_a',line,strg)
3267                          if (strg/="") then
3268                            strg1=trim(strg)
3269                            read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3)
3270                          else
3271                            strg1=trim(line)
3272                            read(strg1,*)(cell(1,1,mu,1,idisp),mu=1,3)
3273                          end if
3274                        end  if
3275                      end if
3276                      if ((line(1:7)==char(60)//'cell_b')) then
3277                        call rdfromline_value('cell_b',line,strg)
3278                        if (strg/="") then
3279                          strg1=trim(strg)
3280                          read(strg1,*) (cell(1,1,mu,2,idisp),mu=1,3)
3281                        else
3282                          read(funit,'(a)',iostat=ios) readline
3283                          call rmtabfromline(readline)
3284                          line=adjustl(readline)
3285                          call rdfromline_value('cell_b',line,strg)
3286                          if (strg/="") then
3287                            strg1=trim(strg)
3288                            read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3)
3289                          else
3290                            strg1=trim(line)
3291                            read(strg1,*)(cell(1,1,mu,2,idisp),mu=1,3)
3292                          end if
3293                        end  if
3294                      end if
3295                      if ((line(1:19)==char(60)//'/displacement_diff')) then
3296                        found2=.false.
3297                      end if
3298                    end do
3299                  end if
3300                end do!end do while displacement
3301 !              In the XML the atom index begin to zero
3302 !              Need to shift for fortran array
3303                atindx(1,1,:,:) = atindx(1,1,:,:) + 1
3304 !              Initialisation of the polynomial_term structure with the values from the
3305 !              previous step
3306                iterm = iterm + 1
3307                call polynomial_term_init(atindx(1,1,:,:),cell(1,1,:,:,:),&
3308 &                                        direction(1,1,:),ndisp,nstrain,terms(1,iterm),&
3309 &                                        power_disp(1,1,:),power_strain(1,1,:),&
3310 &                                        strain(1,1,:),weight(1,1),check=.true.)
3311              end if!end if term
3312            end do!end do while found (coeff)
3313 
3314 !          Initialisation of the polynomial_coefficent structure with the values from the
3315 !          previous step
3316            icoeff = icoeff + 1
3317            call polynomial_coeff_init(coefficient(1),nterm,coeffs(icoeff),terms(1,:))
3318            call polynomial_coeff_getName(name,coeffs(icoeff),symbols,recompute=.true.)
3319            call polynomial_coeff_setName(name,coeffs(icoeff))
3320 !          Deallocation of the terms array for this coefficient
3321            do jj=1,nterm_max
3322              call polynomial_term_free(terms(1,jj))
3323            end do
3324          end if!end if line = coefficient
3325        end if!end if ios==0
3326      end do!end do while on file
3327 
3328      close(unit=funit)
3329 
3330 #endif
3331      ABI_DATATYPE_DEALLOCATE(terms)
3332      ABI_DEALLOCATE(atindx)
3333      ABI_DEALLOCATE(coefficient)
3334      ABI_DEALLOCATE(cell)
3335      ABI_DEALLOCATE(direction)
3336      ABI_DEALLOCATE(strain)
3337      ABI_DEALLOCATE(power_disp)
3338      ABI_DEALLOCATE(power_strain)
3339      ABI_DEALLOCATE(weight)
3340      ABI_DEALLOCATE(symbols)
3341    end if !End if master
3342 
3343 !9-MPI BROADCAST
3344  do ii=1,ncoeff
3345    call polynomial_coeff_broadcast(coeffs(ii),master, comm)
3346  end do
3347 
3348 !10-checks
3349 
3350 !11-debug print
3351  if(debug)then
3352    do ii=1,ncoeff
3353      do jj=1,coeffs(ii)%nterm
3354 #if defined HAVE_LIBXML
3355        write(200+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp
3356        write(200+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx
3357        write(200+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:)
3358        write(200+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:)
3359        write(200+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction
3360        write(200+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp
3361        write(200+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight
3362 #else
3363        write(300+my_rank,*)"ii,jj,ndisp,nterm",ii,jj,coeffs(ii)%nterm,coeffs(ii)%terms(jj)%ndisp
3364        write(300+my_rank,*)"atindx",coeffs(ii)%terms(jj)%atindx
3365        write(300+my_rank,*)"cell1",coeffs(ii)%terms(jj)%cell(:,1,:)
3366        write(300+my_rank,*)"cell2",coeffs(ii)%terms(jj)%cell(:,2,:)
3367        write(300+my_rank,*)"direction",coeffs(ii)%terms(jj)%direction
3368        write(300+my_rank,*)"power_disp",coeffs(ii)%terms(jj)%power_disp
3369        write(300+my_rank,*)"weight",coeffs(ii)%terms(jj)%weight
3370 #endif
3371      end do
3372    end do
3373 #if defined HAVE_LIBXML
3374    close(200+my_rank)
3375 #else
3376    close(300+my_rank)
3377 #endif
3378  end if
3379 
3380 !12-Initialisation of eff_pot
3381  call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff)
3382 
3383 !13-Deallocation of type
3384  do ii=1,ncoeff
3385    call polynomial_coeff_free(coeffs(ii))
3386  end do
3387  ABI_DATATYPE_DEALLOCATE(coeffs)
3388 
3389 end subroutine coeffs_xml2effpot

m_effective_potential_file/effective_potential_file_getDimCoeff [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_getDimCoeff

FUNCTION

 This routine test the xml with polynomial coefficients
 Return the number of coefficients and the maximum number of displacement/strain

INPUTS

 filename = names of the files

OUTPUT

 ncoeff = number of coefficient for the polynome
 nterm(ncoeff) = number terms per coefficient
 ndisp(nterm,ncoeff) = number displacement per term

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

724 subroutine effective_potential_file_getDimCoeff(filename,ncoeff,ndisp_max,nterm_max)
725 
726 
727 !This section has been created automatically by the script Abilint (TD).
728 !Do not modify the following lines by hand.
729 #undef ABI_FUNC
730 #define ABI_FUNC 'effective_potential_file_getDimCoeff'
731 !End of the abilint section
732 
733  implicit none
734 
735 !Arguments ------------------------------------
736 !scalars
737  character(len=fnlen),intent(in) :: filename
738  integer,intent(out) :: ncoeff,ndisp_max,nterm_max
739 !Local variables-------------------------------
740  !scalar
741  integer ::  filetype
742 #ifndef HAVE_LIBXML
743  integer ::  count,count2
744  integer :: funit = 1,ios=0
745  logical :: found,found2
746 #endif
747 !arrays
748 #ifndef HAVE_LIBXML
749  character (len=XML_RECL) :: line,readline
750 #endif
751  character(len=500) :: message
752 
753 ! *************************************************************************
754 
755  call effective_potential_file_getType(filename,filetype)
756 
757  if (filetype==3 .or. filetype==23) then
758    write(message, '(2a)' )' Extraction of the number of coefficient in the XML ',&
759 &                         trim(filename)
760    call wrtout(std_out,message,'COLL')
761 
762    ncoeff = 0
763    nterm_max = 0
764    ndisp_max = 0
765 
766 #if defined HAVE_LIBXML
767 !  Read with libxml the number of coefficient
768    call effpot_xml_getDimCoeff(char_f2c(trim(filename)),ncoeff,nterm_max,ndisp_max)
769 #else
770 !  Read by hand
771 !  Start a reading loop
772    found=.false.
773    ncoeff = 0
774 
775    if (open_file(filename,message,unit=funit,form="formatted",status="old",&
776 &                action="read") /= 0) then
777      MSG_ERROR(message)
778    end if
779 
780 !  First parse to know the number of coefficients
781    ios = 0
782    do while (ios == 0)
783      read(funit,'(a)',iostat=ios) readline
784      if(ios == 0)then
785        call rmtabfromline(readline)
786        line=adjustl(readline)
787 !      Need test with char(9) because the old version of XML file
788 !      from old script includes tarbulation at the begining of each line
789        if (line(1:12)==char(60)//'coefficient') then
790          ncoeff=ncoeff+1
791          count = 0
792          found = .false.
793          do while(.not.found)
794            read(funit,'(a)',iostat=ios) readline
795            call rmtabfromline(readline)
796            line=adjustl(readline)
797            if (line(1:5)==char(60)//'term') then
798              count = count +1
799              found2 = .false.
800              count2 = 0
801              do while(.not.found2)
802                read(funit,'(a)',iostat=ios) readline
803                call rmtabfromline(readline)
804                line=adjustl(readline)
805                if (line(1:13)==char(60)//'displacement') then
806                  count2 = count2 + 1
807                else if (line(1:7)==char(60)//'strain') then
808                  count2 = count2 + 1
809                else if (line(1:6)==char(60)//'/term') then
810                  if (count2 > ndisp_max) ndisp_max = count2
811                  found2 = .true.
812                else
813                  cycle
814                end if
815              end do
816            else  if (line(1:13)==char(60)//'/coefficient') then
817              if (count > nterm_max) nterm_max = count
818              found = .true.
819            else
820              cycle
821            end if
822          end do
823          cycle
824        end if
825      end if
826    end do
827 
828    close(funit)
829 #endif
830 
831  else
832 !  Maybe one day add an other type of file...
833    write(message, '(a,a,a,a)' )&
834 &   ' The file ',trim(filename),' is not compatible with multibinit',ch10
835    MSG_ERROR(message)
836  end if
837 
838 ! Do some checks
839  if (ncoeff < 1) then
840    write(message, '(5a)' )&
841 &   ' Unable to read the number of coeff from ',trim(filename),ch10,&
842 &   ' This file is not compatible with multibinit',ch10
843    MSG_ERROR(message)
844  end if
845 
846 end subroutine effective_potential_file_getDimCoeff

m_effective_potential_file/effective_potential_file_getDimMD [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_getDimMD

FUNCTION

 Read MD FILE (HIST or ASCII) and return the dimensions
 (natom and nstep)

INPUTS

 filename = path of the file

OUTPUT

 natom = number of atoms
 nstep = number of MD steps

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

 981 subroutine effective_potential_file_getDimMD(filename,natom,nstep)
 982 
 983 
 984 !This section has been created automatically by the script Abilint (TD).
 985 !Do not modify the following lines by hand.
 986 #undef ABI_FUNC
 987 #define ABI_FUNC 'effective_potential_file_getDimMD'
 988 !End of the abilint section
 989 
 990  implicit none
 991 
 992 !Arguments ------------------------------------
 993 !scalars
 994  integer,intent(out) :: natom,nstep
 995 !arrays
 996  character(len=fnlen),intent(in) :: filename
 997 !Local variables-------------------------------
 998 !scalar
 999  integer :: ia,natm_old,natm_new
1000  integer :: nenergy,nrprimd
1001  integer :: ios=0,ios2=0,ios3=0
1002  integer :: unit_md=24
1003  logical :: compatible,netcdf
1004 #if defined HAVE_NETCDF
1005  integer :: natom_id,time_id,xyz_id,six_id
1006  integer :: ncid,ncerr
1007  character(len=5) :: char_tmp
1008 #endif
1009 !arrays
1010  character (len=10000) :: readline,line
1011  character(len=500) :: msg
1012 
1013 ! *************************************************************************
1014 
1015  natom = 0
1016  nstep = 0
1017 !try to read netcdf
1018  netcdf = .false.
1019 #if defined HAVE_NETCDF
1020  ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid)
1021  if(ncerr == NF90_NOERR) then
1022    netcdf = .TRUE.
1023    ncerr = nf90_inq_dimid(ncid,"natom",natom_id)
1024    if(ncerr /= NF90_NOERR)  netcdf = .FALSE.
1025    ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id)
1026    if(ncerr /= NF90_NOERR)  netcdf = .FALSE.
1027    ncerr = nf90_inq_dimid(ncid,"time",time_id)
1028    if(ncerr /= NF90_NOERR)  netcdf = .FALSE.
1029    ncerr = nf90_inq_dimid(ncid,"six",six_id)
1030    if(ncerr /= NF90_NOERR)  netcdf = .FALSE.
1031    if(netcdf)then
1032      ncerr = nf90_inquire_dimension(ncid,natom_id,char_tmp,natom)
1033      NCF_CHECK_MSG(ncerr," inquire dimension ID for natom")
1034      ncerr = nf90_inquire_dimension(ncid,time_id,char_tmp,nstep)
1035      NCF_CHECK_MSG(ncerr," inquire dimension ID for time")
1036    end if
1037  end if
1038 #endif
1039 
1040  if(.not.netcdf) then
1041 !  try to read ASCII file...
1042    if (open_file(filename,msg,unit=unit_md,form="formatted",&
1043 &       status="old",action="read") /= 0) then
1044      MSG_ERROR(msg)
1045    end if
1046 
1047 !  Start a reading loop in fortran to get the dimension of the file
1048    rewind(unit=unit_md)
1049    ios = 0
1050    nstep   = 0
1051    nrprimd = 0
1052    natm_old= 0
1053    natm_new= 0
1054    nenergy = 0
1055    compatible = .TRUE.
1056 
1057    do while ((ios==0))
1058 !    special treatment of the first step
1059      if(nstep==0)then
1060        ios2 = 0
1061        do while ((ios2==0))
1062          read(unit_md,'(a)',iostat=ios) readline
1063          line=adjustl(readline)
1064          call elementfromline(line,ia)
1065          if (ia==1)then
1066            nstep = nstep + 1
1067            ios2 = 1
1068          end if
1069        end do
1070      end if
1071      read(unit_md,'(a)',iostat=ios) readline
1072      if(ios == 0)then
1073        line=adjustl(readline)
1074        call elementfromline(line,ia)
1075        if (ia==1)then
1076          nenergy = nenergy + 1
1077          nrprimd = 0
1078        else if(ia==3)then
1079          nrprimd = nrprimd + 1
1080        end if
1081        if(nrprimd == 3)then
1082          ios3 = 0
1083          natm_new = 0
1084          do while ((ios3==0))
1085            read(unit_md,'(a)',iostat=ios3) readline
1086            if(ios3==0)then
1087              line=adjustl(readline)
1088              call elementfromline(line,ia)
1089              if(ia==1)then
1090                if(nstep==1) then
1091                  natm_old = natm_new
1092                else
1093                  if(natm_old /= natm_new) compatible = .FALSE.
1094                end if
1095                ios3 = 1
1096                ios2 = 1
1097                nstep = nstep + 1
1098              end if
1099              if(ia==6)then
1100                natm_new = natm_new + 1
1101              end if
1102            end if!end if ios3
1103          end do
1104        end if ! end if nrprimd
1105      end if! end if os1
1106    end do
1107 
1108    natom = natm_new - 1
1109    if(nstep /= nenergy) compatible = .FALSE.
1110    if(natom <= 0) compatible = .FALSE.
1111    if(nstep <= 0) compatible = .FALSE.
1112 
1113    if(.not.compatible)then
1114      natom = 0
1115      nstep = 0
1116    end if
1117  end if! end if not netcdf
1118 
1119 end subroutine effective_potential_file_getDimMD

m_effective_potential_file/effective_potential_file_getDimStrainCoupling [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_getDimStrainCoupling

FUNCTION

 Return the number of nrpt for specific strain coupling from xml system file

INPUTS

 filename = names of the files
 voigt    = voigt notation of the strain

OUTPUT

 nrpt  = number of rpt points

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

872 subroutine effective_potential_file_getDimStrainCoupling(filename,nrpt,voigt)
873 
874 
875 !This section has been created automatically by the script Abilint (TD).
876 !Do not modify the following lines by hand.
877 #undef ABI_FUNC
878 #define ABI_FUNC 'effective_potential_file_getDimStrainCoupling'
879 !End of the abilint section
880 
881  implicit none
882 
883 !Arguments ------------------------------------
884 !scalars
885  character(len=fnlen),intent(in) :: filename
886  integer,intent(in) :: voigt
887  integer,intent(out) :: nrpt
888 !Local variables-------------------------------
889  !scalar
890 #ifndef HAVE_LIBXML
891  integer :: irpt,ivoigt
892  integer :: funit = 1,ios=0
893  logical :: found
894 #endif
895 !arrays
896 #ifndef HAVE_LIBXML
897  character (len=XML_RECL) :: line,readline,strg,strg1
898  character(len=500) :: message
899 #endif
900 
901 ! *************************************************************************
902 
903    nrpt = 0
904 
905 #if defined HAVE_LIBXML
906 !  Read with libxml the number of coefficient
907    call effpot_xml_getDimStrainCoupling(char_f2c(trim(filename)),nrpt,voigt)
908 #else
909 !  Read by hand
910 !  Start a reading loop
911    found=.false.
912 
913    if (open_file(filename,message,unit=funit,form="formatted",status="old",&
914 &                action="read") /= 0) then
915      MSG_ERROR(message)
916    end if
917 
918 !  First parse to know the number of atoms
919    do while (ios == 0.and.(.not.found))
920      read(funit,'(a)',iostat=ios) readline
921      if(ios == 0)then
922        call rmtabfromline(readline)
923        line=adjustl(readline)
924        if ((line(1:16)=='<strain_coupling')) then
925          read(funit,'(a)',iostat=ios) readline
926          call rdfromline("voigt",line,strg)
927          strg1=trim(strg)
928          read(strg1,*) ivoigt
929          if (ivoigt == voigt)then
930            irpt = 0
931            do while (.not.found)
932              read(funit,'(a)',iostat=ios) readline
933              call rmtabfromline(readline)
934              line=adjustl(readline)
935              if ((line(1:26)=='<correction_force_constant')) then
936                irpt = irpt + 1
937                cycle
938              end if
939              if ((line(1:17)=='</strain_coupling')) then
940                found = .TRUE.
941                nrpt = irpt
942                cycle
943              end if
944            end do
945          else
946            cycle
947          end if
948        end if
949      end if
950    end do
951 
952    close(funit)
953 #endif
954 
955 end subroutine effective_potential_file_getDimStrainCoupling

m_effective_potential_file/effective_potential_file_getDimSystem [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_getDimSystem

FUNCTION

 This routine test the xml or ddb file
 Return the number of atoms/ntypat in the unit cell from ddb and xml
 Return natom/ntypat/nqpt and nrpt if the file is XML file
 In case of DDB file, you have to run bigbx9 to get nrpt

INPUTS

 filename = names of the files

OUTPUT

 natom = number of atoms
 ntypat= number of type of atoms
 nqpt  = number of q points
 nrpt  = number of rpt points

PARENTS

      m_effective_potential_file,multibinit

CHILDREN

SOURCE

597 subroutine effective_potential_file_getDimSystem(filename,natom,ntypat,nqpt,nrpt)
598 
599  use m_ddb
600  use m_ddb_hdr
601 
602 !This section has been created automatically by the script Abilint (TD).
603 !Do not modify the following lines by hand.
604 #undef ABI_FUNC
605 #define ABI_FUNC 'effective_potential_file_getDimSystem'
606 !End of the abilint section
607 
608  implicit none
609 
610 !Arguments ------------------------------------
611 !scalars
612  character(len=fnlen),intent(in) :: filename
613  integer,intent(out) :: natom,ntypat,nqpt,nrpt
614 !arrays
615 
616 !Local variables-------------------------------
617  !scalar
618  integer :: filetype
619 ! integer :: dimekb,lmnmax,mband,mtyp,msym,nblok,nkpt,usepaw
620  integer :: ddbun = 666
621  character(len=500) :: message
622  type(ddb_hdr_type) :: ddb_hdr
623 !arrays
624 ! *************************************************************************
625 
626  natom = 0
627  ntypat= 0
628  nqpt = 0
629  nrpt  = 0
630 
631  call effective_potential_file_getType(filename,filetype)
632 
633  if(filetype==1)then
634    write(message, '(6a)' )ch10,' The file ',trim(filename),ch10,&
635 &                  ' is DDB file (extraction of the number of atoms)',ch10
636    call wrtout(std_out,message,'COLL')
637 
638    write(message, '(8a)' )&
639 &   ' The file ',trim(filename),ch10,' is ddb file only the number of atoms is read,',&
640 &    'if you want to predic the number of cell (nrpt)',ch10,' use bigbx9 routines',ch10
641    call wrtout(std_out,message,'COLL')
642 
643    call ddb_hdr_open_read(ddb_hdr,filename,ddbun,DDB_VERSION,&
644 &                         dimonly=1)
645    natom = ddb_hdr%natom
646    ntypat = ddb_hdr%ntypat
647 
648    call ddb_hdr_free(ddb_hdr)
649 
650 !  Must read some value to initialze  array (nprt for ifc)
651 !   call bigbx9(inp%brav,dummy_cell,0,1,inp%ngqpt,inp%nqshft,nrpt,ddb%rprim,dummy_rpt)
652 
653  else if (filetype==2 .or. filetype==23) then
654    write(message, '(5a)' )ch10,' The file ',trim(filename),&
655 &                ' is XML file (extraction of all informations)'
656    call wrtout(std_out,message,'COLL')
657 
658    call system_getDimFromXML(filename,natom,ntypat,nqpt,nrpt)
659 
660  else
661    write(message, '(a,a,a,a)' )&
662 &   ' The file ',trim(filename),' is not compatible with multibinit',ch10
663    MSG_ERROR(message)
664  end if
665 
666 ! TODO hexu: temporarily disabled. Discuss with alex how to do this properly.
667 ! Do some checks
668 ! if (natom < 1) then
669 !   write(message, '(a,a,a,a,a)' )&
670 !&   ' Unable to read the number of atom from ',trim(filename),ch10,&
671 !&   'This file  is not compatible with multibinit',ch10
672 !   MSG_ERROR(message)
673 ! end if
674 !
675 ! if (filetype==2 .or. filetype==23) then
676 !
677 !   if (natom < 1) then
678 !     write(message, '(a,a,a)' )&
679 !&     ' Unable to read the number of atom from ',trim(filename),ch10
680 !     MSG_ERROR(message)
681 !   end if
682 !
683 !   if (nrpt < 1) then
684 !     write(message, '(a,a,a)' )&
685 !&     ' Unable to read the number of rpt points ',trim(filename),ch10
686 !     MSG_ERROR(message)
687 !   end if
688 !
689 !   if (ntypat < 1) then
690 !     write(message, '(a,a,a)' )&
691 !&     ' Unable to read the number of type of atoms ',trim(filename),ch10
692 !     MSG_ERROR(message)
693 !   end if
694 !
695 ! end if
696 
697 end subroutine effective_potential_file_getDimSystem

m_effective_potential_file/effective_potential_file_getType [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_getType

FUNCTION

 This routine test the xml or ddb file

INPUTS

 filename = names of the files

OUTPUT

 type_file  = 0 no type found
              1 DDB file
              2 XML file the system definition and harmonic part
              3 XML file with polynomial coefficients
             23 XML file with both system definition and polynomial coefficients
             40 NetCDF file with history of MD or snapshot
             41 ASCII file with history of MD or snapshot

PARENTS

      m_effective_potential_file,multibinit

CHILDREN

SOURCE

464 subroutine effective_potential_file_getType(filename,filetype)
465 
466 
467 !This section has been created automatically by the script Abilint (TD).
468 !Do not modify the following lines by hand.
469 #undef ABI_FUNC
470 #define ABI_FUNC 'effective_potential_file_getType'
471 !End of the abilint section
472 
473  implicit none
474 
475 !Arguments ------------------------------------
476 !scalars
477  character(len=fnlen),intent(in) :: filename
478  integer, intent(out) :: filetype
479 !arrays
480 !Local variables-------------------------------
481 !scalar
482  integer :: natom,nstep
483  integer :: ddbun = 666,ios=0
484  character(len=500) :: message
485  character (len=1000) :: line,readline
486 #if defined HAVE_NETCDF
487  integer :: natom_id,time_id,xyz_id,six_id
488  integer :: ncid,ncerr
489  logical :: md_file
490 #endif
491 
492 !arrays
493 ! *************************************************************************
494 
495  filetype = 0
496 
497  if (open_file(filename,message,unit=ddbun,form="formatted",status="old",action="read") /= 0) then
498    MSG_ERROR(message)
499  end if
500 
501 !Check if the file is a XML file or a DDB and in this case, store the DDB code.
502  ios = 0
503  do while ((ios==0))
504    read(ddbun,'(a)',iostat=ios) readline
505    call rmtabfromline(readline)
506    line=adjustl(readline)
507    if(line(3:13)=="xml version") then
508      do while ((ios==0))
509        read(ddbun,'(a)',iostat=ios) readline
510        call rmtabfromline(readline)
511        line=adjustl(readline)
512        if(line(1:16)==char(60)//"Heff_definition".or.&
513 &         line(1:17)==char(60)//"Terms_definition")then
514          filetype = 3
515          ios = -1
516        end if
517        if(line(1:18)==char(60)//"System_definition") then
518          filetype = 2
519          do while ((ios==0))
520            read(ddbun,'(a)',iostat=ios) readline
521            call rmtabfromline(readline)
522            line=adjustl(readline)
523            if(line(1:16)==char(60)//"Heff_definition".or.&
524 &             line(1:17)==char(60)//"Terms_definition")then
525              filetype = 23
526              ios = -1
527            end if
528          end do
529        end if
530      end do
531    else  if(line(6:24)=="DERIVATIVE DATABASE") then
532      filetype = 1
533      ios = -1
534    end if
535  end do
536  close(ddbun)
537 
538  if(filetype/=0) return
539 
540 !try to read netcdf HIST file
541 #if defined HAVE_NETCDF
542  ncerr=nf90_open(path=trim(filename),mode=NF90_NOWRITE,ncid=ncid)
543  if(ncerr == NF90_NOERR) then
544    md_file = .TRUE.
545    ncerr = nf90_inq_dimid(ncid,"natom",natom_id)
546    if(ncerr /= NF90_NOERR)  md_file = .FALSE.
547    ncerr = nf90_inq_dimid(ncid,"xyz",xyz_id)
548    if(ncerr /= NF90_NOERR)  md_file = .FALSE.
549    ncerr = nf90_inq_dimid(ncid,"time",time_id)
550    if(ncerr /= NF90_NOERR)  md_file = .FALSE.
551    ncerr = nf90_inq_dimid(ncid,"six",six_id)
552    if(ncerr /= NF90_NOERR)  md_file = .FALSE.
553    if (md_file) then
554      filetype = 40
555      return
556    end if
557  end if
558  ncerr = nf90_close(ncid)
559 #endif
560 
561  if(filetype/=0) return
562 
563 !Try to get the dim of MD ASCII file
564  call effective_potential_file_getDimMD(filename,natom,nstep)
565  if(natom /= 0 .and. nstep/=0) filetype = 41
566 
567 end subroutine effective_potential_file_getType

m_effective_potential_file/effective_potential_file_mapHistToRef [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_mapHistToRef

FUNCTION

 Generate the supercell in the effective potential according to the size of the
 supercell in the hist file
 Check if the hist file match to reference supercell in the effective potential
 If not, the hist file is reordering

INPUTS

 eff_pot<type(effective_potential)> = effective potential
 hist<type(abihist)> = The history of the MD
 comm = MPI communicator

OUTPUT

 hist<type(abihist)> = The history of the MD

PARENTS

      m_fit_polynomial_coeff,multibinit

CHILDREN

SOURCE

3539 subroutine effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose)
3540 
3541 
3542 !This section has been created automatically by the script Abilint (TD).
3543 !Do not modify the following lines by hand.
3544 #undef ABI_FUNC
3545 #define ABI_FUNC 'effective_potential_file_mapHistToRef'
3546 !End of the abilint section
3547 
3548  implicit none
3549 
3550 !Arguments ------------------------------------
3551 !scalars
3552  integer,intent(in) :: comm
3553  logical,optional,intent(in) :: verbose
3554 !arrays
3555  type(effective_potential_type),intent(inout) :: eff_pot
3556  type(abihist),intent(inout) :: hist
3557 !Local variables-------------------------------
3558 !scalar
3559  integer :: factE_hist,ia,ib,ii,jj,natom_hist,ncells,nstep_hist
3560  real(dp):: factor
3561  logical :: revelant_factor,need_map,need_verbose
3562 !arrays
3563  real(dp) :: rprimd_hist(3,3),rprimd_ref(3,3),scale_cell(3)
3564  integer :: ncell(3)
3565  integer,allocatable  :: blkval(:),list(:)
3566  real(dp),allocatable :: xred_hist(:,:),xred_ref(:,:)
3567  character(len=500) :: msg
3568  type(abihist) :: hist_tmp
3569 ! *************************************************************************
3570 
3571 !Set optional values
3572  need_verbose = .false.
3573  if (present(verbose)) need_verbose = verbose
3574 
3575  natom_hist = size(hist%xred,2)
3576  nstep_hist = size(hist%xred,3)
3577 
3578 !Try to set the supercell according to the hist file
3579  rprimd_ref(:,:)  = eff_pot%crystal%rprimd
3580  rprimd_hist(:,:) = hist%rprimd(:,:,1)
3581 
3582  do ia=1,3
3583    scale_cell(:) = zero
3584    do ii=1,3
3585      if(abs(rprimd_ref(ii,ia)) > tol10)then
3586        scale_cell(ii) = rprimd_hist(ii,ia) / rprimd_ref(ii,ia)
3587      end if
3588    end do
3589 !  Check if the factor for the supercell is revelant
3590    revelant_factor = .TRUE.
3591    do ii=1,3
3592      if(abs(scale_cell(ii)) < tol10) cycle
3593      factor = abs(scale_cell(ii))
3594      do jj=ii,3
3595        if(abs(scale_cell(jj)) < tol10) cycle
3596        if(abs(abs(scale_cell(ii))-abs(scale_cell(jj))) > tol10) revelant_factor = .FALSE.
3597      end do
3598    end do
3599    if(.not.revelant_factor)then
3600      write(msg, '(3a)' )&
3601 &         'unable to map the hist file ',ch10,&
3602 &         'Action: check/change your MD file'
3603      MSG_ERROR(msg)
3604    else
3605      ncell(ia) = nint(factor)
3606    end if
3607  end do
3608 
3609  ncells = product(ncell)
3610 
3611 !Check if the energy store in the hist is revelant, sometimes some MD files gives
3612 !the energy of the unit cell... This is not suppose to happen... But just in case...
3613  do ii=1,nstep_hist
3614    factE_hist = int(anint(hist%etot(ii) / eff_pot%energy))
3615    if(factE_hist == 1) then
3616 !    In this case we mutiply the energy of the hist by the number of cell
3617      hist%etot(ii) = hist%etot(ii)  * ncells
3618    end if
3619    if(factE_hist /=1 .and. factE_hist /= ncells)then
3620      write(msg, '(4a,I0,a,I0,2a,I0,3a,I0,3a)' )ch10,&
3621 &          ' --- !WARNING',ch10,&
3622 &          '     The energy of the step ',ii,' seems to be with multiplicity of ',factE_hist,ch10,&
3623 &          '     However, the multiplicity of the cell is ',ncells,'.',ch10,&
3624 &          '     Please check the energy of the step ',ii,ch10,&
3625 &          ' ---',ch10
3626      if(need_verbose) call wrtout(std_out,msg,'COLL')
3627    end if
3628  end do
3629 
3630 
3631 !Set the new supercell datatype into the effective potential reference
3632  call effective_potential_setSupercell(eff_pot,comm,ncell)
3633 
3634 !allocation
3635  ABI_ALLOCATE(blkval,(natom_hist))
3636  ABI_ALLOCATE(list,(natom_hist))
3637  ABI_ALLOCATE(xred_hist,(3,natom_hist))
3638  ABI_ALLOCATE(xred_ref,(3,natom_hist))
3639  blkval = 1
3640  list   = 0
3641  call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3642 &                eff_pot%supercell%xcart,xred_ref)
3643 
3644  xred_hist = hist%xred(:,:,1)
3645 
3646  if(need_verbose) then
3647    write(msg,'(2a,I2,a,I2,a,I2)') ch10,&
3648 &       ' The size of the supercell for the fit is ',ncell(1),' ',ncell(2),' ',ncell(3)
3649    call wrtout(std_out,msg,'COLL')
3650    call wrtout(ab_out,msg,'COLL')
3651  end if
3652 
3653 !try to map
3654  do ia=1,natom_hist
3655    do ib=1,natom_hist
3656      if(blkval(ib)==1)then
3657        if(abs((xred_ref(1,ia)-xred_hist(1,ib))) < 0.1 .and.&
3658 &         abs((xred_ref(2,ia)-xred_hist(2,ib))) < 0.1 .and.&
3659 &         abs((xred_ref(3,ia)-xred_hist(3,ib))) < 0.1) then
3660          blkval(ib) = 0
3661          list(ib) = ia
3662        end if
3663      end if
3664    end do
3665  end do
3666 
3667 !Check before transfert
3668  if(.not.all(blkval==0))then
3669    write(msg, '(5a)' )&
3670 &         'Unable to map the molecular dynamic file ',ch10,&
3671 &         'on the reference supercell structure',ch10,&
3672 &         'Action: change the MD file'
3673      MSG_ERROR(msg)
3674  end if
3675 
3676  do ia=1,natom_hist
3677    if(.not.any(list(:)==ia))then
3678      write(msg, '(5a)' )&
3679 &         'Unable to map the molecular dynamic file  ',ch10,&
3680 &         'on the reference supercell structure',ch10,&
3681 &         'Action: change the MD file'
3682      MSG_ERROR(msg)
3683    end if
3684  end do
3685 
3686  need_map = .FALSE.
3687  do ia=1,natom_hist
3688    if(list(ia) /= ia) need_map = .TRUE.
3689  end do
3690  if(need_map)then
3691    if(need_verbose) then
3692      write(msg, '(11a)' )ch10,&
3693 &      ' --- !WARNING',ch10,&
3694 &      '     The ordering of the atoms in the hist file is different,',ch10,&
3695 &      '     of the one built by multibinit. The hist file will be map,',ch10,&
3696 &      '     on the ordering of multibinit.',ch10,&
3697 &      ' ---',ch10
3698      call wrtout(ab_out,msg,'COLL')
3699      call wrtout(std_out,msg,'COLL')
3700    end if
3701 
3702 ! Allocate hist datatype
3703    call abihist_init(hist_tmp,natom_hist,nstep_hist,.false.,.false.)
3704 ! copy all the information
3705    do ia=1,nstep_hist
3706      hist%ihist = ia
3707      hist_tmp%ihist = ia
3708      call abihist_copy(hist,hist_tmp)
3709    end do
3710    hist_tmp%mxhist = nstep_hist
3711 
3712 ! reoder array
3713    do ia=1,natom_hist
3714      hist_tmp%xred(:,list(ia),:)=hist%xred(:,ia,:)
3715      hist_tmp%fcart(:,list(ia),:)=hist%fcart(:,ia,:)
3716      hist_tmp%vel(:,list(ia),:)=hist%vel(:,ia,:)
3717    end do
3718 
3719 ! free the old hist and reinit
3720    call abihist_free(hist)
3721    call abihist_init(hist,natom_hist,nstep_hist,.false.,.false.)
3722 ! copy the temporary hist into output
3723    do ia=1,nstep_hist
3724      hist%ihist = ia
3725      hist_tmp%ihist = ia
3726      call abihist_copy(hist_tmp,hist)
3727    end do
3728    hist_tmp%mxhist = nstep_hist
3729 
3730    call abihist_free(hist_tmp)
3731  end if
3732 
3733 !deallocation
3734  ABI_DEALLOCATE(blkval)
3735  ABI_DEALLOCATE(list)
3736  ABI_DEALLOCATE(xred_hist)
3737  ABI_DEALLOCATE(xred_ref)
3738 
3739 end subroutine effective_potential_file_mapHistToRef

m_effective_potential_file/effective_potential_file_readMDfile [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 effective_potential_file_readMDfile

FUNCTION

 Read MD FILE (HIST or ASCII)

INPUTS

 filename = path of the file
 option,optional   = 0 (default), the stress is printed in the MD File
                     1, the force on the cell is printed in the MD File (-1 * stress),
                        in this case, we multiply the stress by -1 in order to get the stresse

OUTPUT

 hist<type(abihist)> = datatype with the  history of the MD

PARENTS

      m_effective_potential_file,multibinit

CHILDREN

SOURCE

3416 subroutine effective_potential_file_readMDfile(filename,hist,option)
3417 
3418 
3419 !This section has been created automatically by the script Abilint (TD).
3420 !Do not modify the following lines by hand.
3421 #undef ABI_FUNC
3422 #define ABI_FUNC 'effective_potential_file_readMDfile'
3423 !End of the abilint section
3424 
3425  implicit none
3426 
3427 !Arguments ------------------------------------
3428 !scalars
3429  integer,optional :: option
3430 !arrays
3431  type(abihist),intent(inout) :: hist
3432  character(len=fnlen),intent(in) :: filename
3433 !Local variables-------------------------------
3434 !scalar
3435  integer :: ia,ii,mu,nu,natom,nstep,type,option_in
3436  integer :: ios=0, unit_md=24
3437 !arrays
3438  character (len=10000) :: readline,line
3439  real(dp) :: tmp(6)
3440  real(dp),allocatable :: xcart(:,:)
3441 
3442 ! *************************************************************************
3443 
3444  call effective_potential_file_getType(filename,type)
3445 
3446  option_in = 0
3447  if(present(option))then
3448    option_in = option
3449  end if
3450 
3451  if(type==40)then
3452 !  Netcdf type
3453    call read_md_hist(filename,hist,.FALSE.,.FALSE.,.FALSE.)
3454 
3455  else if(type==41)then
3456 
3457 !  ASCII file
3458    call  effective_potential_file_getDimMD(filename,natom,nstep)
3459 
3460    ii  = 1
3461    ios = 0
3462 
3463    ABI_ALLOCATE(xcart,(3,natom))
3464    call abihist_free(hist)
3465    call abihist_init(hist,natom,nstep,.FALSE.,.FALSE.)
3466 
3467 !  Start a reading loop in fortran
3468    rewind(unit=unit_md)
3469    do while ((ios==0).and.ii<=nstep)
3470      read(unit_md,'(a)',iostat=ios) readline
3471      read(unit_md,'(a)',iostat=ios) readline
3472      line=adjustl(readline)
3473      read(line,*) hist%etot(ii)
3474      hist%etot(ii) = hist%etot(ii)
3475      do mu=1,3
3476        read(unit_md,'(a)',iostat=ios) readline
3477        line=adjustl(readline)
3478        read(line,*) (hist%rprimd(nu,mu,ii),nu=1,3)
3479      end do
3480      do ia=1,natom
3481        read(unit_md,'(a)',iostat=ios) readline
3482        line=adjustl(readline)
3483        read(line,*) (tmp(mu),mu=1,6)
3484        xcart(:,ia) = tmp(1:3)
3485        hist%fcart(:,ia,ii) = tmp(4:6)
3486      end do
3487      call xcart2xred(natom,hist%rprimd(:,:,ii),xcart(:,:),hist%xred(:,:,ii))
3488      read(unit_md,'(a)',iostat=ios) readline
3489      line=adjustl(readline)
3490      read(line,*) (hist%strten(mu,ii),mu=1,6)
3491      ii = ii + 1
3492    end do
3493    do ii=1,nstep
3494      do mu=1,3
3495        hist%acell(mu,:) = hist%rprimd(mu,mu,ii)
3496      end do
3497    end do
3498    close(unit_md)
3499    ABI_DEALLOCATE(xcart)
3500 
3501  end if!end if type
3502 
3503    if((type==40 .or. type==41).and.option == 1)then
3504 !    multiply by -1 if the current strten -1*stress, we need only stress...
3505      hist%strten(:,:) = -1 * hist%strten(:,:)
3506    end if
3507 
3508 
3509 
3510 end subroutine effective_potential_file_readMDfile

m_effective_potential_file/elementfromline [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 elementfromline

FUNCTION

 Read the number of element of a line

INPUTS

  line= string from which the data are read

OUTPUT

  nelement = number of element in the line

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

3836 subroutine elementfromline(line,nelement)
3837 
3838 
3839 !This section has been created automatically by the script Abilint (TD).
3840 !Do not modify the following lines by hand.
3841 #undef ABI_FUNC
3842 #define ABI_FUNC 'elementfromline'
3843 !End of the abilint section
3844 
3845  implicit none
3846 
3847 !Arguments ---------------------------------------------
3848  character(len=*), intent(in) :: line
3849  integer, intent(out) :: nelement
3850 !Local variables ---------------------------------------
3851  integer :: ii,n
3852  logical :: element
3853 ! *********************************************************************
3854 
3855 !Set the output
3856  nelement = 0
3857  n = len_trim(line)
3858  element = .false.
3859  do ii=1,n
3860    if(.not.element.and.line(ii:ii) /="")  then
3861      element=.true.
3862    else
3863      if((element.and.line(ii:ii) =="")) then
3864        element=.false.
3865        nelement = nelement + 1
3866      end if
3867    end if
3868    if((element.and.ii==n)) nelement = nelement + 1
3869  end do
3870 
3871  end subroutine elementfromline

m_effective_potential_file/rdfromline [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 rdfromline

FUNCTION

 Read the value of a keyword from a XML line
 Same function than m_pawxmlps/paw_rdfromline.F90

INPUTS

  keyword= keyword which value has to be read
  line= string from which the data are read (line from a XML)

OUTPUT

  output= (string) value of the keyword

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

3896  subroutine rdfromline(keyword,line,output)
3897 
3898 
3899 !This section has been created automatically by the script Abilint (TD).
3900 !Do not modify the following lines by hand.
3901 #undef ABI_FUNC
3902 #define ABI_FUNC 'rdfromline'
3903 !End of the abilint section
3904 
3905  implicit none
3906 
3907 !Arguments ---------------------------------------------
3908   character(len=*), intent(in) :: keyword,line
3909   character(len=*), intent(out) :: output
3910 !Local variables ---------------------------------------
3911   character(len=len(line)) :: temp
3912   integer :: pos,pos2
3913 
3914 ! *********************************************************************
3915 
3916  output=""
3917  pos=index(line,trim(keyword))
3918  if (pos>0) then
3919    temp=line(pos+len_trim(keyword):len_trim(line))
3920    pos=index(temp,char(34))
3921    if (pos>0) then
3922      pos2=index(temp(pos+1:len_trim(temp)),char(34))
3923      if (pos2>0) then
3924        output=temp(pos+1:pos+pos2-1)
3925      end if
3926    end if
3927  end if
3928 
3929  end subroutine rdfromline

m_effective_potential_file/rdfromline_value [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 rdfromline

FUNCTION

 Read the value of a keyword from a XML line

INPUTS

  keyword= keyword which value has to be read
  line= string from which the data are read (line from a XML)

OUTPUT

  output= (string) value of the keyword

PARENTS

      system_xml2effpot

CHILDREN

      paw_rdfromline

SOURCE

4005  subroutine rdfromline_value(keyword,line,output)
4006 
4007 
4008 !This section has been created automatically by the script Abilint (TD).
4009 !Do not modify the following lines by hand.
4010 #undef ABI_FUNC
4011 #define ABI_FUNC 'rdfromline_value'
4012 !End of the abilint section
4013 
4014  implicit none
4015 
4016 !Arguments ---------------------------------------------
4017   character(len=*), intent(in) :: keyword,line
4018   character(len=*), intent(out) :: output
4019 !Local variables ---------------------------------------
4020   character(len=len(line)) :: temp
4021   integer :: pos,pos2
4022 
4023 ! *********************************************************************
4024 
4025  output=""
4026  pos=index(line,trim(keyword))
4027  if (pos==2) then
4028    pos=pos+len_trim(keyword)
4029    pos=pos+index(line(pos:len_trim(line)),char(62))
4030    temp=line(pos:len_trim(line))
4031    pos2=index(temp,char(60))
4032    if (pos2>0) then
4033      output=line(pos:pos+pos2-2)
4034    else
4035      output=line(pos:len_trim(line))
4036    end if
4037  else
4038    if(pos>2)then
4039      output=line(1:pos-3)
4040    end if
4041  end if
4042 end subroutine rdfromline_value

m_effective_potential_file/rmtabfromline [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 rmtabfromline

FUNCTION

 Read remove tab from the begining of line

INPUTS

  line= string from which the data are read (line from a XML)

OUTPUT

  output= line without tab

PARENTS

    system_xml2effpot

CHILDREN

    rmtabfromline

SOURCE

3954 recursive subroutine rmtabfromline(line)
3955 
3956 
3957 !This section has been created automatically by the script Abilint (TD).
3958 !Do not modify the following lines by hand.
3959 #undef ABI_FUNC
3960 #define ABI_FUNC 'rmtabfromline'
3961 !End of the abilint section
3962 
3963  implicit none
3964 
3965 !Arguments ---------------------------------------------
3966   character(len=*), intent(inout) :: line
3967 !Local variables ---------------------------------------
3968   integer :: pos
3969 
3970 ! *********************************************************************
3971 
3972  pos=index(line,char(9))
3973  if (pos==1) then
3974    line = line(2:len_trim(line))//" "
3975    call rmtabfromline(line)
3976  end if
3977 
3978  end subroutine rmtabfromline

m_effective_potential_file/system_ddb2effpot [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 system_ddb2effpot

FUNCTION

  Transfert ddb into effective potential structure.
  Also calculate the IFC

INPUTS

 crytal<type(crystal_t)> = datatype with all the information for the crystal
 ddb<type(ddb_type)> = datatype with the ddb
 inp<type(multibinit_dtset_type)> = datatype with the input variables of multibinit
 comm = MPI communicator

OUTPUT

 effective_potantial<type(effective_potential_type)> = effective_potential datatype to be initialized

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

2285 subroutine system_ddb2effpot(crystal,ddb, effective_potential,inp,comm)
2286 
2287  use defs_basis
2288  use m_errors
2289  use m_abicore
2290  use m_dynmat
2291  use m_xmpi
2292 
2293  use m_ddb
2294  use m_ifc
2295  use m_copy,            only : alloc_copy
2296  use m_crystal,         only : crystal_t,crystal_print
2297  use m_multibinit_dataset, only : multibinit_dtset_type
2298  use m_effective_potential, only : effective_potential_type, effective_potential_free
2299 
2300 !This section has been created automatically by the script Abilint (TD).
2301 !Do not modify the following lines by hand.
2302 #undef ABI_FUNC
2303 #define ABI_FUNC 'system_ddb2effpot'
2304 !End of the abilint section
2305 
2306  implicit none
2307 
2308 !Arguments ------------------------------------
2309 !scalars
2310  integer,intent(in) :: comm
2311 !arrays
2312  type(ddb_type),intent(inout) :: ddb
2313  type(effective_potential_type), intent(inout) :: effective_potential
2314  type(crystal_t),intent(in) :: crystal
2315  type(multibinit_dtset_type),intent(in) :: inp
2316 
2317 !Local variables-------------------------------
2318 !scalar
2319  real(dp):: wcount1,wcount2
2320  integer :: chneut,i1,i2,i3,ia,ib,iblok,idir1,idir2,ierr,ii,ipert1,iphl1
2321  integer :: ipert2,irpt,irpt2,ivarA,ivarB,max1,max2,max3,min1,min2,min3
2322  integer :: msize,mpert,natom,nblok,nrpt_new,nrpt_new2,rftyp,selectz
2323  integer :: my_rank,nproc,prt_internalstr
2324  logical :: iam_master=.FALSE.
2325  integer,parameter :: master=0
2326  integer :: nptsym,nsym
2327  integer :: msym = 192,  use_inversion = 1, space_group
2328  real(dp):: max_phfq,tolsym = tol8
2329 !arrays
2330  integer :: bravais(11),cell_number(3),cell2(3)
2331  integer :: shift(3),rfelfd(4),rfphon(4),rfstrs(4)
2332  integer,allocatable :: cell_red(:,:)
2333  real(dp):: dielt(3,3),elast_clamped(6,6),fact
2334  real(dp):: red(3,3),qphnrm(3),qphon(3,3)
2335  real(dp),allocatable :: blkval(:,:,:,:,:,:),d2asr(:,:,:,:,:)
2336  real(dp),allocatable :: instrain(:,:),zeff(:,:,:)
2337  real(dp),pointer :: atmfrc_red(:,:,:,:,:),wghatm_red(:,:,:)
2338  character(len=500) :: message
2339  type(asrq0_t) :: asrq0
2340  type(ifc_type) :: ifc
2341  real(dp),allocatable :: d2cart(:,:,:,:,:),displ(:)
2342  real(dp),allocatable :: eigval(:,:),eigvec(:,:,:,:,:),phfrq(:)
2343  real(dp),allocatable :: spinat(:,:),tnons(:,:)
2344  integer,allocatable  :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:)
2345 
2346 ! *************************************************************************
2347 
2348 !0 MPI variables
2349  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2350  iam_master = (my_rank == master)
2351 
2352 !Free the eff_pot before filling
2353  call effective_potential_free(effective_potential)
2354 
2355 !Initialisation of usefull values
2356   natom = ddb%natom
2357   nblok = ddb%nblok
2358   mpert=natom+6
2359   msize=3*mpert*3*mpert;
2360 
2361 !Tranfert the ddb into usable array (ipert and idir format like in abinit)
2362   ABI_ALLOCATE(blkval,(2,3,mpert,3,mpert,nblok))
2363   blkval = 0
2364   blkval = reshape(ddb%val,(/2,3,mpert,3,mpert,nblok/))
2365 
2366 !**********************************************************************
2367 ! Transfert crystal values
2368 !**********************************************************************
2369 ! Re-generate symmetry operations from the lattice and atomic coordinates
2370   ABI_ALLOCATE(spinat,(3,natom))
2371   ABI_ALLOCATE(ptsymrel,(3,3,msym))
2372   ABI_ALLOCATE(symafm,(msym))
2373   ABI_ALLOCATE(symrel,(3,3,msym))
2374   ABI_ALLOCATE(tnons,(3,msym))
2375   spinat = zero;  symrel = 0;  symafm = 0;  tnons = zero ; space_group = 0;
2376   call symlatt(bravais,msym,nptsym,ptsymrel,crystal%rprimd,tolsym)
2377   call symfind(0,(/zero,zero,zero/),crystal%gprimd,0,msym,crystal%natom,0,nptsym,nsym,&
2378 &              0,0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,&
2379 &              crystal%typat,use_inversion,crystal%xred)
2380   if(crystal%nsym/=nsym)then
2381     write(message,'(4a,I0,3a,I0,3a)') ch10,&
2382 &          ' --- !WARNING:',ch10,&
2383 &          '     There is ',nsym,' found for the crystal',ch10,&
2384 &          '     but ',crystal%nsym,' found in the DDB',ch10,&
2385 &          ' ---'
2386       call wrtout(std_out,message,'COLL')
2387   end if
2388 
2389   call crystal_init(ddb%amu,effective_potential%crystal,&
2390 &                   space_group,crystal%natom,crystal%npsp,&
2391 &                   crystal%ntypat,nsym,crystal%rprimd,&
2392 &                   crystal%typat,crystal%xred,crystal%zion,&
2393 &                   crystal%znucl,crystal%timrev,crystal%use_antiferro,&
2394 &                   .FALSE.,crystal%title,&
2395 &                   symrel=symrel,tnons=tnons,&
2396 &                   symafm=symafm)
2397 
2398   ABI_DEALLOCATE(spinat)
2399   ABI_DEALLOCATE(ptsymrel)
2400   ABI_DEALLOCATE(symafm)
2401   ABI_DEALLOCATE(symrel)
2402   ABI_DEALLOCATE(tnons)
2403 
2404 !**********************************************************************
2405 ! Transfert energy from input file
2406 !**********************************************************************
2407   write(message, '(2a,(80a),6a)') ch10,('=',ii=1,80),ch10,ch10,&
2408 &     ' Extraction of the energy of the structure (unit: Hartree)',ch10
2409   call wrtout(std_out,message,'COLL')
2410   call wrtout(ab_out,message,'COLL')
2411   if (ddb_get_etotal(ddb,effective_potential%energy) == 0) then
2412     if(abs(inp%energy_reference) < tol16)then
2413       write(message,'(5a)')&
2414 &      ' Warning : Energy of the reference structure is not specify in',&
2415 &      ' the input file.',ch10,' Energy will set to zero',ch10
2416       call wrtout(std_out,message,'COLL')
2417       effective_potential%energy = zero
2418     else
2419       effective_potential%energy = inp%energy_reference
2420     end if
2421   else
2422     if(abs(inp%energy_reference) > tol16)then
2423       write(message,'(6a)')&
2424 &      ' Warning : Energy of the reference structure is specify in',&
2425 &      ' the input file.',ch10,' and in the DDB.',&
2426 &      ' The value of the energy is set with the value from the input file',ch10
2427       call wrtout(std_out,message,'COLL')
2428       effective_potential%energy = inp%energy_reference
2429     end if
2430   end if
2431   write(message,'(a,es25.12)') ' Energy  = ',&
2432 &                    effective_potential%energy
2433   call wrtout(std_out,message,'COLL')
2434   call wrtout(ab_out,message,'COLL')
2435 
2436 !**********************************************************************
2437 ! Dielectric Tensor and Effective Charges
2438 !**********************************************************************
2439   ABI_ALLOCATE(zeff,(3,3,natom))
2440   ABI_ALLOCATE(effective_potential%harmonics_terms%zeff,(3,3,natom))
2441 
2442   rftyp   = 1 ! Blocks obtained by a non-stationary formulation.
2443   chneut  = 1 ! The ASR for effective charges is imposed
2444   selectz = 0 ! No selection of some parts of the effective charge tensor
2445   iblok = ddb_get_dielt_zeff(ddb,crystal,rftyp,chneut,selectz,dielt,zeff)
2446   if (iblok /=0 .and. maxval(abs(dielt)) < 10000) then
2447     effective_potential%harmonics_terms%epsilon_inf = dielt
2448     effective_potential%harmonics_terms%zeff = zeff
2449   else
2450     effective_potential%harmonics_terms%epsilon_inf(1,1) = one
2451     effective_potential%harmonics_terms%epsilon_inf(2,2) = one
2452     effective_potential%harmonics_terms%epsilon_inf(3,3) = one
2453     effective_potential%harmonics_terms%zeff = zero
2454   end if
2455 
2456 !**********************************************************************
2457 ! Look after the blok no. that contains the stress tensor
2458 !**********************************************************************
2459   write(message, '(a,a,(80a),a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2460 &   ' Extraction of the stress tensor (unit: GPa) and forces (unit: Ha/bohr)'
2461   call wrtout(std_out,message,'COLL')
2462   call wrtout(ab_out,message,'COLL')
2463 
2464   ABI_ALLOCATE(effective_potential%fcart,(3,natom))
2465   effective_potential%fcart = zero
2466   effective_potential%strten = zero
2467 
2468   qphon(:,1)=zero
2469   qphnrm(1)=zero
2470   rfphon(1:2)=0
2471   rfelfd(1:2)=0
2472   rfstrs(1:2)=0
2473   rftyp=4
2474 
2475   call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2476 
2477   if (iblok /=0) then
2478    if(any(abs(inp%strten_reference)>tol16))then
2479      write(message,'(10a)') ch10,&
2480 &          ' --- !WARNING:',ch10,&
2481 &          '     The stress tensor of the reference structure is specify in the',ch10,&
2482 &          '     input file and in the DDB. The value of the stress tensor is set',ch10,&
2483 &          '     with the value from the input file',ch10,&
2484 &          ' ---'
2485      call wrtout(std_out,message,'COLL')
2486      call wrtout(ab_out,message,'COLL')
2487      effective_potential%strten = inp%strten_reference
2488    else
2489 !    firts give the corect stress values store in hartree
2490 !    diagonal parts
2491      effective_potential%strten(1)=blkval(1,1,natom+3,1,1,iblok) *  crystal%ucvol
2492      effective_potential%strten(2)=blkval(1,2,natom+3,1,1,iblok) *  crystal%ucvol
2493      effective_potential%strten(3)=blkval(1,3,natom+3,1,1,iblok) *  crystal%ucvol
2494 !    the shear parts
2495      effective_potential%strten(4)=blkval(1,1,natom+4,1,1,iblok) *  crystal%ucvol
2496      effective_potential%strten(5)=blkval(1,2,natom+4,1,1,iblok) *  crystal%ucvol
2497      effective_potential%strten(6)=blkval(1,3,natom+4,1,1,iblok) *  crystal%ucvol
2498    end if
2499 !  Get forces
2500    effective_potential%fcart(:,1:natom) = blkval(1,:,1:natom,1,1,iblok)
2501  else
2502    if(all(abs(inp%strten_reference(:))<tol16))then
2503      write(message,'(8a)') ch10,&
2504 &          ' --- !WARNING:',ch10,&
2505 &          '     The stress tensor of the reference structure is not specify',ch10,&
2506 &          '     The stress tensor will be set to zero',ch10,&
2507 &          ' ---'
2508      call wrtout(std_out,message,'COLL')
2509      call wrtout(ab_out,message,'COLL')
2510      effective_potential%strten = zero
2511    else
2512      effective_potential%strten = inp%strten_reference
2513    end if
2514  end if
2515 
2516  if(any(abs(effective_potential%strten(:)) >tol16))then
2517    write(message, '(3a)' )ch10,&
2518 &   ' Cartesian components of forces (hartree/bohr)',ch10
2519    call wrtout(ab_out,message,'COLL')
2520    call wrtout(std_out,  message,'COLL')
2521    do ii = 1, natom
2522      write(message, '(I4,a,3(e16.8))' ) &
2523 &     ii,'   ',effective_potential%fcart(:,ii)
2524 
2525      call wrtout(ab_out,message,'COLL')
2526      call wrtout(std_out,  message,'COLL')
2527    end do
2528 
2529    write(message, '(a,a)' )ch10,&
2530 &   ' Cartesian components of stress tensor (hartree/bohr^3)'
2531    call wrtout(ab_out,message,'COLL')
2532    call wrtout(std_out,  message,'COLL')
2533    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2534 &   '  sigma(1 1)=',effective_potential%strten(1) / crystal%ucvol,&
2535 &   '  sigma(3 2)=',effective_potential%strten(4) / crystal%ucvol
2536    call wrtout(ab_out,message,'COLL')
2537    call wrtout(std_out,  message,'COLL')
2538    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2539 &   '  sigma(2 2)=',effective_potential%strten(2) / crystal%ucvol,&
2540 &   '  sigma(3 1)=',effective_potential%strten(5) / crystal%ucvol
2541    call wrtout(ab_out,message,'COLL')
2542    call wrtout(std_out,  message,'COLL')
2543    write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2544 &   '  sigma(3 3)=',effective_potential%strten(3) / crystal%ucvol,&
2545 &   '  sigma(2 1)=',effective_potential%strten(6) / crystal%ucvol
2546    call wrtout(ab_out,message,'COLL')
2547    call wrtout(std_out,  message,'COLL')
2548    write(message, '(a)' ) ' '
2549    call wrtout(ab_out,message,'COLL')
2550    call wrtout(std_out,  message,'COLL')
2551  end if
2552 
2553 !**********************************************************************
2554 ! Elastic tensors at Gamma Point
2555 !**********************************************************************
2556   write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2557 &   ' Extraction of the clamped elastic tensor (unit:10^2GPa)',ch10
2558   call wrtout(std_out,message,'COLL')
2559   call wrtout(ab_out,message,'COLL')
2560 
2561 ! look after the blok no.iblok that contains the elastic tensor
2562   qphon(:,1)=zero
2563   qphnrm(1)=zero
2564   rfphon(1:2)=0
2565   rfelfd(1:2)=0
2566   rfstrs(1:2)=3 ! Need uniaxial  both stresses and  shear stresses
2567   rftyp=1 ! Blocks obtained by a non-stationary formulation.
2568 ! for both diagonal and shear parts
2569   call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2570 
2571   if (iblok /=0) then
2572 !   extraction of the elastic constants from the blkvals (GPa)
2573     do ivarA=1,6
2574       do ivarB=1,6
2575 !       because the elastic constant is 6*6,
2576 !       so we should judge if the idir is larger than 3
2577 !       or not
2578         if(ivarA>3) then
2579           idir1=ivarA-3
2580           ipert1=natom+4  !for the shear modulus
2581         else if(ivarA<=3) then
2582           idir1=ivarA
2583           ipert1=natom+3  !for the diagonal part
2584         end if
2585         if(ivarB>3) then
2586           idir2=ivarB-3
2587           ipert2=natom+4  !for the shear modulus
2588         else if(ivarB<=3) then
2589           idir2=ivarB
2590           ipert2=natom+3  !for the diagonal part
2591         end if
2592         elast_clamped(ivarA,ivarB) = blkval(1,idir1,ipert1,idir2,ipert2,iblok)
2593       end do
2594     end do
2595     fact=HaBohr3_GPa / crystal%ucvol
2596     do ivarA=1,6
2597       write(message,'(6f12.7)')elast_clamped(ivarA,1)*fact/100.00_dp,&
2598 &                              elast_clamped(ivarA,2)*fact/100.00_dp,&
2599 &                              elast_clamped(ivarA,3)*fact/100.00_dp,&
2600 &                              elast_clamped(ivarA,4)*fact/100.00_dp,&
2601 &                              elast_clamped(ivarA,5)*fact/100.00_dp,&
2602 &                              elast_clamped(ivarA,6)*fact/100.00_dp
2603     call wrtout(std_out,message,'COLL')
2604     call wrtout(ab_out,message,'COLL')
2605     end do
2606 
2607 !   Set the clamped tensor into the effective potentiel
2608     effective_potential%harmonics_terms%elastic_constants = elast_clamped
2609 
2610   else
2611 
2612     write(message,'(3a)')ch10,&
2613 &    ' Warning : Elastic Tensor is set to zero (not available in the DDB)'
2614     call wrtout(std_out,message,'COLL')
2615     call wrtout(ab_out,message,'COLL')
2616 
2617 !   Set the clamped tensor to zero into the effective potentiel (not available in the DDB)
2618     effective_potential%harmonics_terms%elastic_constants = zero
2619   end if
2620 
2621 !**********************************************************************
2622 !   Acoustic Sum Rule
2623 !***************************************************************************
2624 ! ASR-correction (d2asr) has to be determined here from the Dynamical matrix at Gamma.
2625   ABI_ALLOCATE(d2asr,(2,3,natom,3,natom))
2626 
2627   write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2628 &   ' Calculation of acoustic sum rule',ch10
2629   call wrtout(std_out,message,'COLL')
2630   call wrtout(ab_out,message,'COLL')
2631 
2632 ! Find the Gamma block in the DDB (no need for E-field entries)
2633   qphon(:,1)=zero
2634   qphnrm(1)=zero
2635   rfphon(1:2)=1
2636   rfelfd(:)=0
2637   rfstrs(:)=0
2638   rftyp=inp%rfmeth
2639 
2640   call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2641 
2642   d2asr = zero
2643   if (iblok /=0) then
2644     call asria_calc(inp%asr,d2asr,ddb%val(:,:,iblok),ddb%mpert,ddb%natom)
2645   end if
2646 
2647   ! Acoustic Sum Rule
2648   ! In case the interatomic forces are not calculated, the
2649   ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma.
2650   asrq0 = ddb_get_asrq0(ddb, inp%asr, inp%rfmeth, crystal%xcart)
2651 
2652 !**********************************************************************
2653 ! Interatomic Forces Calculation
2654 !**********************************************************************
2655 ! ifc to be calculated for interpolation
2656   write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,&
2657 &   ' Calculation of the interatomic forces from DDB',ch10
2658   call wrtout(std_out,message,'COLL')
2659   call wrtout(ab_out,message,'COLL')
2660 
2661   call ifc_init(ifc,crystal,ddb,inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,&
2662 &   inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,effective_potential%harmonics_terms%zeff,&
2663 &   inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm)
2664 
2665 !***************************************************************************
2666 ! Interpolation of the dynamical matrix for each qpoint from ifc
2667 !***************************************************************************
2668 
2669   ABI_ALLOCATE(d2cart,(2,3,mpert,3,mpert))
2670   ABI_ALLOCATE(displ,(2*3*natom*3*natom))
2671   ABI_ALLOCATE(eigval,(3,natom))
2672   ABI_ALLOCATE(eigvec,(2,3,natom,3,natom))
2673   ABI_ALLOCATE(phfrq,(3*natom))
2674 
2675   ABI_ALLOCATE(effective_potential%harmonics_terms%dynmat,(2,3,natom,3,natom,inp%nph1l))
2676   ABI_ALLOCATE(effective_potential%harmonics_terms%phfrq,(3*natom,inp%nph1l))
2677   ABI_ALLOCATE(effective_potential%harmonics_terms%qpoints,(3,inp%nph1l))
2678 
2679   write(message,'(a,(80a),3a)')ch10,('=',ii=1,80),ch10,ch10,&
2680 &     ' Calculation of dynamical matrix for each ph1l points '
2681   call wrtout(ab_out,message,'COLL')
2682   call wrtout(std_out,message,'COLL')
2683 
2684 !Transfer value in effective_potential structure
2685   effective_potential%harmonics_terms%nqpt         = inp%nph1l
2686   effective_potential%harmonics_terms%qpoints(:,:) = inp%qph1l(:,:)
2687 
2688 ! Store the highest frequency
2689   max_phfq = zero
2690 
2691   do iphl1=1,inp%nph1l
2692 
2693    ! Initialisation of the phonon wavevector
2694     qphon(:,1)=inp%qph1l(:,iphl1)
2695     if (inp%nph1l /= 0) qphnrm(1) = inp%qnrml1(iphl1)
2696 
2697     ! Get d2cart using the interatomic forces and the
2698     ! long-range coulomb interaction through Ewald summation
2699     call gtdyn9(ddb%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip,ifc%dyewq0,d2cart,crystal%gmet,&
2700 &     ddb%gprim,mpert,natom,ifc%nrpt,qphnrm(1),qphon(:,1),crystal%rmet,ddb%rprim,ifc%rpt,&
2701 &     ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,zeff)
2702 
2703     ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
2704     call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,crystal%indsym,&
2705 &     mpert,crystal%nsym,natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon,&
2706 &     crystal%rprimd,inp%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol)
2707 
2708     ! Write the phonon frequencies
2709     call dfpt_prtph(displ,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
2710 
2711 !   Store the highest frequency in cmm-1
2712     max_phfq = max(maxval(phfrq*Ha_cmm1),max_phfq)
2713 
2714     effective_potential%harmonics_terms%dynmat(:,:,:,:,:,iphl1) = d2cart(:,:,:natom,:,:natom)
2715     effective_potential%harmonics_terms%phfrq(:,iphl1) = phfrq(:) * Ha_cmm1
2716 
2717   end do
2718 
2719   write(message, '(2a,f15.7,a)' ) ch10,&
2720 &   ' The highest frequency found is ',max_phfq,' cm-1'
2721   call wrtout(std_out,message,'COLL')
2722 
2723   ABI_DEALLOCATE(d2cart)
2724   ABI_DEALLOCATE(displ)
2725   ABI_DEALLOCATE(eigval)
2726   ABI_DEALLOCATE(eigvec)
2727   ABI_DEALLOCATE(phfrq)
2728 
2729 !**********************************************************************
2730 ! Transfert inter-atomic forces constants in reduced coordinates
2731 !**********************************************************************
2732 
2733 !Reorder cell from canonical coordinates to reduced coordinates (for multibinit)
2734 !store the number of ifc before rearrangement
2735 
2736 ! Store the sum of the weight of IFC for the final check
2737   wcount1 = 0
2738   do irpt=1,ifc%nrpt
2739     wcount1 = wcount1 + sum(ifc%wghatm(:,:,irpt))
2740   end do
2741 
2742 !Set the maximum and the miminum for the bound of the cell
2743   max1 = maxval(ifc%cell(1,:));  min1 = minval(ifc%cell(1,:))
2744   max2 = maxval(ifc%cell(2,:));  min2 = minval(ifc%cell(2,:))
2745   max3 = maxval(ifc%cell(3,:));  min3 = minval(ifc%cell(3,:))
2746   cell_number(1) = max1 - min1 + 1
2747   cell_number(2) = max2 - min2 + 1
2748   cell_number(3) = max3 - min3 + 1
2749 
2750 ! set the new number of cell, sometimes, in canonical coordinates,
2751 ! some cell are delete but they exist in reduced coordinates.
2752   nrpt_new = product(cell_number(:))
2753 
2754 ! Allocate temporary array
2755   ABI_ALLOCATE(atmfrc_red,(3,natom,3,natom,nrpt_new))
2756   ABI_ALLOCATE(wghatm_red,(natom,natom,nrpt_new))
2757   ABI_ALLOCATE(cell_red,(3,nrpt_new))
2758 
2759   wghatm_red(:,:,:) = zero
2760 
2761   if(iam_master)then
2762     do ia=1,natom
2763       do ib=1,natom
2764 
2765 !       Simple Lattice
2766         if (inp%brav==1) then
2767 !          In this case, it is better to work in reduced coordinates
2768 !          As rcan is in canonical coordinates, => multiplication by gprim
2769            do ii=1,3
2770              red(1,ii)=  ifc%rcan(1,ia)*ddb%gprim(1,ii) + &
2771  &                       ifc%rcan(2,ia)*ddb%gprim(2,ii) + &
2772  &                       ifc%rcan(3,ia)*ddb%gprim(3,ii)
2773              red(2,ii)=  ifc%rcan(1,ib)*ddb%gprim(1,ii) + &
2774                          ifc%rcan(2,ib)*ddb%gprim(2,ii) + &
2775  &                       ifc%rcan(3,ib)*ddb%gprim(3,ii)
2776            end do
2777          end if
2778 
2779 !       Get the shift of cell
2780         shift(:) = int(anint(red(2,:) - crystal%xred(:,ib)) - anint(red(1,:) - crystal%xred(:,ia)))
2781 
2782         do irpt=1,ifc%nrpt
2783 
2784           cell2(:)= int(ifc%cell(:,irpt) + shift(:))
2785 
2786 !         Use boundary condition to get the right cell
2787           if (cell2(1) < min1 .and. cell2(1) < max1) then
2788             cell2(1) = cell2(1) + cell_number(1)
2789           else if (cell2(1) > min1 .and. cell2(1) > max1) then
2790             cell2(1) = cell2(1) - cell_number(1)
2791           end if
2792 
2793           if (cell2(2) < min2 .and. cell2(2) < max2) then
2794             cell2(2) = cell2(2) + cell_number(2)
2795           else if (cell2(2) > min2 .and. cell2(2) > max2) then
2796             cell2(2) = cell2(2) - cell_number(2)
2797           end if
2798 
2799           if (cell2(3) < min3 .and. cell2(3) < max3) then
2800             cell2(3) = cell2(3) + cell_number(3)
2801           else if (cell2(3) > min3 .and. cell2(3) > max3) then
2802             cell2(3) = cell2(3) - cell_number(3)
2803           end if
2804 
2805            irpt2=1
2806            do i1=min1,max1
2807              do i2=min2,max2
2808                do i3=min3,max3
2809                  if (i1  ==  cell2(1)  .and.&
2810                      i2  ==  cell2(2)  .and.&
2811                      i3  ==  cell2(3)) then
2812                    wghatm_red(ia,ib,irpt2) =  ifc%wghatm(ia,ib,irpt)
2813                    atmfrc_red(:,ia,:,ib,irpt2) = ifc%atmfrc(:,ia,:,ib,irpt)
2814                    cell_red(1,irpt2) = i1
2815                    cell_red(2,irpt2) = i2
2816                    cell_red(3,irpt2) = i3
2817                  end if
2818                  irpt2 = irpt2 + 1
2819                end do
2820              end do
2821           end do
2822         end do
2823       end do
2824     end do
2825   end if
2826 
2827   call xmpi_bcast(atmfrc_red,master, comm, ierr)
2828   call xmpi_bcast(wghatm_red,master, comm, ierr)
2829   call xmpi_bcast(cell_red,master, comm, ierr)
2830 
2831 !  Copy ifc into effective potential
2832 ! !!Warning eff_pot%ifcs only contains atmfrc,short_atmfrc,ewald_atmfrc,,nrpt and cell!!
2833 ! rcan,ifc%rpt,wghatm and other quantities
2834 ! are not needed for effective potential!!!
2835   call ifc_free(ifc)
2836   call ifc_free(effective_potential%harmonics_terms%ifcs)
2837 
2838 ! Only conserve the necessary points in rpt
2839   nrpt_new2 = 0
2840   do irpt = 1, nrpt_new
2841     if (abs(sum(wghatm_red(:,:,irpt))) >tol16) then
2842       nrpt_new2 = nrpt_new2 + 1
2843     end if
2844   end do
2845 
2846 ! Set the new number of rpt
2847   effective_potential%harmonics_terms%ifcs%nrpt = nrpt_new2
2848 
2849 ! Allocation of the final arrays
2850   ABI_ALLOCATE(effective_potential%harmonics_terms%ifcs%atmfrc,(3,natom,3,natom,nrpt_new2))
2851   ABI_ALLOCATE(effective_potential%harmonics_terms%ifcs%short_atmfrc,(3,natom,3,natom,nrpt_new2))
2852   ABI_ALLOCATE(effective_potential%harmonics_terms%ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt_new2))
2853   ABI_ALLOCATE(effective_potential%harmonics_terms%ifcs%cell,(3,nrpt_new2))
2854   ABI_ALLOCATE(effective_potential%harmonics_terms%ifcs%wghatm,(natom,natom,nrpt_new2))
2855 
2856   irpt2 = 0
2857   do irpt = 1,nrpt_new
2858     if (abs(sum(wghatm_red(:,:,irpt))) > tol16) then
2859       irpt2 = irpt2 + 1
2860 !     Apply weight on each R point
2861       do ia=1,effective_potential%crystal%natom
2862         do ib=1,effective_potential%crystal%natom
2863           atmfrc_red(:,ia,:,ib,irpt) = atmfrc_red(:,ia,:,ib,irpt)*wghatm_red(ia,ib,irpt)
2864         end do
2865       end do
2866       effective_potential%harmonics_terms%ifcs%cell(:,irpt2) = cell_red(:,irpt)
2867       effective_potential%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt2) = atmfrc_red(:,:,:,:,irpt)
2868       if (inp%dipdip == 1) then
2869         effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=&
2870 &                                                                     atmfrc_red(:,:,:,:,irpt)
2871       else
2872         effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2) = zero
2873       end if
2874       effective_potential%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt2)=atmfrc_red(:,:,:,:,irpt)
2875       effective_potential%harmonics_terms%ifcs%ewald_atmfrc(:,:,:,:,irpt2) = zero
2876       effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt2) =  wghatm_red(:,:,irpt)
2877     end if
2878   end do
2879 
2880 
2881   ABI_DEALLOCATE(atmfrc_red)
2882   ABI_DEALLOCATE(wghatm_red)
2883   ABI_DEALLOCATE(cell_red)
2884 
2885 ! Final check
2886   wcount2 = 0
2887   do irpt = 1, effective_potential%harmonics_terms%ifcs%nrpt
2888     wcount2 = wcount2 + sum(effective_potential%harmonics_terms%ifcs%wghatm(:,:,irpt))
2889   end do
2890 
2891   if (abs(wcount1-wcount2)/(wcount1+wcount2)>tol8) then
2892     write(message,'(2a,es15.4,a,es15.4,a,es15.4)')'The total wghatm has changed',ch10,&
2893 &    wcount1,' before and ', wcount2, ' now, difference being ',wcount1-wcount2
2894     MSG_BUG(message)
2895   end if
2896 
2897 
2898 !**********************************************************************
2899 ! Internal strain tensors at Gamma point
2900 !**********************************************************************
2901   write(message, '(a,a,(80a),a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
2902 &   ' Calculation of the internal-strain  tensor'
2903   call wrtout(std_out,message,'COLL')
2904   call wrtout(ab_out,message,'COLL')
2905   ABI_ALLOCATE(instrain,(3*natom,6))
2906 ! looking after the no. of blok that contains the internal strain tensor
2907   qphon(:,1)=zero
2908   qphnrm(1)=zero
2909   rfphon(1:2)=0
2910   rfelfd(1:2)=0
2911   rfstrs(1:2)=3
2912   rftyp=1
2913   call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp)
2914 
2915   ABI_ALLOCATE(effective_potential%harmonics_terms%strain_coupling,(6,3,natom))
2916   effective_potential%harmonics_terms%strain_coupling = zero
2917 
2918   if (iblok /=0) then
2919 
2920 !   then print the internal strain tensor (only the force one)
2921     prt_internalstr=1
2922     call ddb_internalstr(inp%asr,ddb%val,d2asr,iblok,instrain,&
2923 &                        ab_out,mpert,natom,nblok,prt_internalstr)
2924 
2925     do ipert1=1,6
2926       do ipert2=1,natom
2927         do idir2=1,3
2928           ii=3*(ipert2-1)+idir2
2929             effective_potential%harmonics_terms%strain_coupling(ipert1,idir2,ipert2)=&
2930 &                                                            (-1.0_dp)*instrain(ii,ipert1)
2931         end do
2932       end do
2933     end do
2934   else
2935     write(message,'(3a)')ch10,&
2936 &    ' Warning : Internal strain is set to zero (not available in the DDB)'
2937     call wrtout(std_out,message,'COLL')
2938     call wrtout(ab_out,message,'COLL')
2939   end if
2940 !-------------------------------------------------------------------------------------
2941 ! DEALLOCATION OF ARRAYS
2942   ABI_DEALLOCATE(blkval)
2943   ABI_DEALLOCATE(zeff)
2944   ABI_DEALLOCATE(instrain)
2945   ABI_DEALLOCATE(d2asr)
2946   call asrq0_free(asrq0)
2947 
2948   write(message,'(a)')ch10
2949   call wrtout(std_out,message,'COLL')
2950   call wrtout(ab_out,message,'COLL')
2951 
2952 end subroutine system_ddb2effpot

m_effective_potential_file/system_getDimFromXML [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 system_getDimFromXML

FUNCTION

 Open xml file of effective potentiel, then reads the variables that
 must be known in order to dimension the arrays before complete reading

INPUTS

 character(len=*) filnam: name of input or output file

OUTPUT

 natom=number of atoms
 ntypat=number of atom types
 nrpt  =number of real space points used to integrate IFC

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

1146 subroutine system_getDimFromXML(filename,natom,ntypat,nph1l,nrpt)
1147 
1148 
1149 !This section has been created automatically by the script Abilint (TD).
1150 !Do not modify the following lines by hand.
1151 #undef ABI_FUNC
1152 #define ABI_FUNC 'system_getDimFromXML'
1153 !End of the abilint section
1154 
1155  implicit none
1156 
1157  !Arguments ------------------------------------
1158  !scalars
1159   character(len=fnlen),intent(in) :: filename
1160   integer, intent(out) :: natom,ntypat,nph1l,nrpt
1161  !arrays
1162  !Local variables-------------------------------
1163  !scalar
1164   integer :: nrpt1,nrpt2
1165   real :: itypat
1166   character(len=500) :: message
1167 #ifndef HAVE_LIBXML
1168   integer :: funit = 1,ios = 0
1169   integer :: iatom
1170   logical :: found
1171   character (len=XML_RECL) :: line,readline
1172   character (len=XML_RECL) :: strg,strg1
1173 #endif
1174   !arrays
1175 #ifndef HAVE_LIBXML
1176   integer,allocatable :: typat(:)
1177 #endif
1178 
1179  ! *************************************************************************
1180 
1181 !Open the atomicdata XML file for reading
1182  write(message,'(5a)') ' system_getDimFromXML :',&
1183 &    '-Opening the file ',trim(filename),' to read dimensions',&
1184 &    ' (before initialisation)'
1185 
1186  call wrtout(std_out,message,'COLL')
1187 
1188  natom = 0
1189  ntypat= 0
1190  nph1l = 0
1191  nrpt  = 0
1192  nrpt1 = 0
1193  nrpt2 = 0
1194  itypat= 0
1195 
1196 !Open the atomicdata XML file for reading
1197 
1198 #if defined HAVE_LIBXML
1199 !Read with libxml
1200  call effpot_xml_getDimSystem(char_f2c(trim(filename)),natom,ntypat,nph1l,nrpt1,nrpt2)
1201 #else
1202 !Read by hand
1203 
1204 !Start a reading loop
1205  found=.false.
1206 
1207  if (open_file(filename,message,unit=funit,form="formatted",status="old",&
1208 &              action="read") /= 0) then
1209    MSG_ERROR(message)
1210  end if
1211 
1212 !First parse to know the number of atoms
1213  do while ((ios==0).and.(.not.found))
1214    read(funit,'(a)',iostat=ios) readline
1215    if(ios ==0)then
1216      call rmtabfromline(readline)
1217      line=adjustl(readline)
1218 
1219 !Need test with char(9) because the old version of XML file
1220 !from old script includes tarbulation at the begining of each line
1221      if (line(1:5)==char(60)//'atom') then
1222        natom=natom+1
1223        cycle
1224      end if
1225 
1226      if (line(1:21)==char(60)//'local_force_constant') then
1227        nrpt1 = nrpt1+1
1228        cycle
1229      end if
1230 
1231      if (line(1:21)==char(60)//'total_force_constant') then
1232        nrpt2 = nrpt2+1
1233        cycle
1234      end if
1235 
1236      if (line(1:7)==char(60)//'qpoint') then
1237        nph1l =  nph1l+1
1238        cycle
1239      end if
1240    end if
1241  end do
1242 
1243 !second parse to get the number of typat
1244  ABI_ALLOCATE(typat,(natom))
1245  typat = 0
1246  iatom = 0
1247 
1248  rewind(funit)
1249 !Start a reading loop
1250  ios   = 0
1251  found = .false.
1252 
1253  do while ((ios==0).and.(.not.found))
1254    read(funit,'(a)',iostat=ios) readline
1255    if(ios == 0)then
1256      call rmtabfromline(readline)
1257      line=adjustl(readline)
1258 
1259      if (line(1:5)==char(60)//'atom') then
1260        iatom = iatom + 1
1261        call rdfromline("mass",line,strg)
1262        strg1=trim(strg)
1263        read(unit=strg1,fmt=*) itypat
1264        if (.not.any(typat==int(itypat))) then
1265          ntypat= ntypat+1
1266        end if
1267        typat(iatom) = int(itypat)
1268      end if
1269 
1270      if (line(1:6)==char(60)//'local') then
1271        found=.true.
1272      end if
1273    end if
1274  end do
1275 
1276  close(funit)
1277  ABI_DEALLOCATE(typat)
1278 
1279 #endif
1280 
1281 !Check the RPT
1282  if (nrpt2/=nrpt1) then
1283    if(nrpt1> 0 .and. nrpt2== 0) then
1284      continue;
1285    else if (nrpt1==0.and.nrpt2>=0) then
1286      write(message, '(5a)' )ch10,&
1287 &   ' WARNING: the number of local IFC is set to 0  ',ch10,&
1288 &   '          Dipdip must be set to zero',ch10
1289      call wrtout(std_out,message,'COLL')
1290    else if (nrpt2 > nrpt1) then
1291      write(message, '(2a,I0,3a,I0,5a)' )ch10,&
1292 &   ' WARNING: the number of total IFC  (',nrpt2,') is not equal to the  ',ch10,&
1293 &   '          the number of short range IFC (',nrpt1,') in ',filename,ch10,&
1294 &   '          the missing ifc will be set to zero',ch10
1295      call wrtout(std_out,message,'COLL')
1296    else if(nrpt1>nrpt2)then
1297      write(message, '(2a,I0,3a,I0,5a)' )ch10,&
1298 &   ' The number of total IFC  (',nrpt2,') is inferior to  ',ch10,&
1299 &   ' the number of short range IFC (',nrpt1,') in ',filename,ch10,&
1300 &   ' This is not possible',ch10
1301      MSG_BUG(message)
1302    end if
1303  end if
1304 
1305 !nrpt is the max between local and total:
1306  nrpt = max(nrpt1,nrpt2)
1307 
1308 
1309 end subroutine system_getDimFromXML

m_effective_potential_file/system_xml2effpot [ Functions ]

[ Top ] [ m_effective_potential_file ] [ Functions ]

NAME

 system_xml2effpot

FUNCTION

 Open xml file of effective potentiel, then reads the variables
 and store them in effective potentential type

INPUTS

 eff_pot<type(effective_potential_type)> = datatype with all the informations for effective potential
 comm=MPI communicator
 character(len=*) filnam: name of input or output file
 strcpling = optional,logical to disable the strcpling

OUTPUT

 eff_pot<type(effective_potential_type)> = datatype with all the informations for effective potential

PARENTS

      m_effective_potential_file

CHILDREN

SOURCE

1336  subroutine system_xml2effpot(eff_pot,filename,comm,strcpling)
1337 
1338  use m_atomdata
1339  use m_effective_potential, only : effective_potential_type
1340  use m_multibinit_dataset, only : multibinit_dtset_type
1341  use m_ab7_symmetry
1342 
1343 !This section has been created automatically by the script Abilint (TD).
1344 !Do not modify the following lines by hand.
1345 #undef ABI_FUNC
1346 #define ABI_FUNC 'system_xml2effpot'
1347 !End of the abilint section
1348 
1349  implicit none
1350 
1351  !Arguments ------------------------------------
1352  !scalars
1353  character(len=*),intent(in) :: filename
1354  integer, intent(in) :: comm
1355  integer, optional,intent(in) :: strcpling
1356  !arrays
1357  type(effective_potential_type), intent(inout) :: eff_pot
1358 
1359  !Local variables-------------------------------
1360  !scalar
1361  integer :: ierr,ii,itypat,my_rank,msym,natom,ncoeff,nrpt,nrpt_scoupling
1362  integer :: ntypat,nph1l,nptsym,npsp,nproc,nsym,space_group,timrev,use_inversion,voigt
1363  real(dp):: energy,tolsym,ucvol
1364  character(len=500) :: message
1365  integer,parameter :: master=0
1366  logical :: has_anharmonics = .FALSE.
1367  logical :: iam_master
1368 #ifndef HAVE_LIBXML
1369  integer :: funit = 1,ios=0
1370  integer :: iatom,iamu,iph1l,irpt,irpt1,irpt2,irpt3,jj,mu,nu
1371  real(dp):: amu
1372  logical :: found,found2,short_range,total_range
1373  character (len=XML_RECL) :: line,readline
1374  character (len=XML_RECL) :: strg,strg1
1375  logical :: has_straincoupling = .FALSE.
1376 #endif
1377  !arrays
1378  integer :: bravais(11)
1379  integer,allocatable :: typat(:)
1380  integer,allocatable  :: symrel(:,:,:),symafm(:),ptsymrel(:,:,:)
1381  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1382  real(dp) :: elastic_constants(6,6),elastic3rd(6,6,6),epsilon_inf(3,3)
1383  real(dp),allocatable :: all_amu(:),cell_local(:,:),cell_total(:,:)
1384  real(dp),allocatable :: elastic_displacement(:,:,:,:),dynmat(:,:,:,:,:,:)
1385  real(dp),allocatable :: local_atmfrc(:,:,:,:,:),total_atmfrc(:,:,:,:,:)
1386  real(dp),allocatable :: spinat(:,:),strain_coupling(:,:,:),phfrq(:,:),qph1l(:,:),tnons(:,:)
1387  real(dp),allocatable :: xcart(:,:),xred(:,:),zeff(:,:,:),znucl(:),zion(:)
1388  character(len=132),allocatable :: title(:)
1389  type(ifc_type) :: ifcs
1390  type(ifc_type),dimension(:),allocatable :: phonon_strain
1391  type(crystal_t)  :: crystal
1392  type(atomdata_t) :: atom
1393 #ifdef HAVE_LIBXML
1394  real(dp),allocatable :: phonon_strain_atmfrc(:,:,:,:,:)
1395  integer,allocatable  :: phonon_straincell(:,:)
1396 #endif
1397 #ifndef HAVE_LIBXML
1398  real(dp),allocatable :: work2(:,:)
1399 #endif
1400 
1401 ! *************************************************************************
1402 
1403 
1404  !Open the atomicdata XML file for reading
1405  write(message,'(a,a)')'-Opening the file ',filename
1406 
1407  call wrtout(ab_out,message,'COLL')
1408  call wrtout(std_out,message,'COLL')
1409 
1410  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1411  iam_master = (my_rank == master)
1412 
1413 !Get Dimention of system and allocation/initialisation of array
1414  call effective_potential_file_getDimSystem(filename,natom,ntypat,nph1l,nrpt)
1415  gmet= zero; gprimd = zero; rmet = zero; rprimd = zero
1416  elastic_constants = zero; epsilon_inf = zero; ncoeff = 0
1417  ABI_ALLOCATE(all_amu,(ntypat))
1418  ABI_ALLOCATE(cell_local,(3,nrpt))
1419  ABI_ALLOCATE(cell_total,(3,nrpt))
1420  ABI_ALLOCATE(elastic_displacement,(6,6,3,natom))
1421  ABI_ALLOCATE(ifcs%atmfrc,(3,natom,3,natom,nrpt))
1422  ABI_ALLOCATE(ifcs%cell,(3,nrpt))
1423  ABI_ALLOCATE(ifcs%short_atmfrc,(3,natom,3,natom,nrpt))
1424  ABI_ALLOCATE(ifcs%ewald_atmfrc,(3,natom,3,natom,nrpt))
1425  ABI_ALLOCATE(strain_coupling,(6,3,natom))
1426  ABI_ALLOCATE(total_atmfrc,(3,natom,3,natom,nrpt))
1427  ABI_ALLOCATE(local_atmfrc,(3,natom,3,natom,nrpt))
1428  ABI_ALLOCATE(dynmat,(2,3,natom,3,natom,nph1l))
1429  ABI_ALLOCATE(typat,(natom))
1430  ABI_ALLOCATE(phfrq,(3*natom,nph1l))
1431  ABI_ALLOCATE(qph1l,(3,nph1l))
1432  ABI_ALLOCATE(xcart,(3,natom))
1433  ABI_ALLOCATE(xred,(3,natom))
1434  ABI_ALLOCATE(zeff,(3,3,natom))
1435  ABI_ALLOCATE(zion,(ntypat))
1436  ABI_ALLOCATE(znucl,(ntypat))
1437 
1438  ABI_DATATYPE_ALLOCATE(phonon_strain,(6))
1439  nrpt_scoupling = 0
1440  do ii = 1,6
1441 !  Get The size of the strainPhonon-coupling
1442    call effective_potential_file_getDimStrainCoupling(filename,nrpt_scoupling,ii-1)
1443    ABI_ALLOCATE(phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt_scoupling))
1444    ABI_ALLOCATE(phonon_strain(ii)%cell,(3,nrpt_scoupling))
1445    phonon_strain(ii)%nrpt   = nrpt_scoupling
1446    phonon_strain(ii)%atmfrc = zero
1447    phonon_strain(ii)%cell   = 0
1448  end do
1449 
1450  all_amu(:) = zero
1451  dynmat(:,:,:,:,:,:)  = zero
1452  cell_local(:,:) = 99D99
1453  cell_total(:,:) = 99D99
1454  elastic3rd(:,:,:) = zero
1455  elastic_displacement(:,:,:,:) = zero
1456  ifcs%nrpt = nrpt
1457  ifcs%atmfrc(:,:,:,:,:)  = zero
1458  ifcs%cell(:,:)  = 0
1459  ifcs%ewald_atmfrc(:,:,:,:,:) = zero
1460  ifcs%short_atmfrc(:,:,:,:,:) = zero
1461  strain_coupling(:,:,:) = zero
1462  phfrq = zero
1463  qph1l = 0
1464  xcart = zero
1465  zeff  = zero
1466  znucl = zero
1467 
1468  if(iam_master)then
1469 !Open the atomicdata XML file for reading
1470 #if defined HAVE_LIBXML
1471 
1472    write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),&
1473 &   ' with LibXML library'
1474 
1475    call wrtout(ab_out,message,'COLL')
1476    call wrtout(std_out,message,'COLL')
1477 
1478 !  Read with libxml library
1479    call effpot_xml_readSystem(char_f2c(trim(filename)),natom,ntypat,nrpt,nph1l,all_amu,&
1480 &                       ifcs%atmfrc,ifcs%cell,dynmat,elastic_constants,energy,&
1481 &                       epsilon_inf,ifcs%ewald_atmfrc,phfrq,rprimd,qph1l,&
1482 &                       ifcs%short_atmfrc,typat,xcart,zeff)
1483 
1484 !  convert atomic mass unit to znucl
1485    do itypat=1,ntypat
1486      do ii=1,103
1487        call atomdata_from_znucl(atom,real(ii,dp))
1488        if (abs((real(atom%amu,sp)-real(all_amu(itypat),sp))&
1489 &         /real(all_amu(itypat),sp)*100)<0.1) then
1490          znucl(itypat) = atom%znucl
1491          exit
1492        end if
1493      end do
1494    end do
1495 
1496 !  Get the Phonon Strain coupling
1497    do voigt = 1,6
1498      nrpt_scoupling = phonon_strain(voigt)%nrpt
1499      ABI_ALLOCATE(phonon_straincell,(3,nrpt_scoupling))
1500      ABI_ALLOCATE(phonon_strain_atmfrc,(3,natom,3,natom,nrpt_scoupling))
1501 
1502 !      Get The value
1503        call effpot_xml_readStrainCoupling(char_f2c(trim(filename)),natom,nrpt_scoupling,(voigt-1),&
1504 &                                         elastic3rd(voigt,:,:),elastic_displacement(voigt,:,:,:),&
1505 &                                         strain_coupling(voigt,:,:),&
1506 &                                         phonon_strain_atmfrc,phonon_straincell)
1507 
1508 !      Check if the 3rd order strain_coupling is present
1509        if(any(elastic3rd>tol10).or.any(elastic_displacement>tol10)) has_anharmonics = .TRUE.
1510        phonon_strain(voigt)%atmfrc(:,:,:,:,:) = phonon_strain_atmfrc(:,:,:,:,:)
1511        phonon_strain(voigt)%cell(:,:)   = phonon_straincell(:,:)
1512        if(any(phonon_strain(voigt)%atmfrc > tol10)) has_anharmonics = .TRUE.
1513 
1514        ABI_DEALLOCATE(phonon_straincell)
1515        ABI_DEALLOCATE(phonon_strain_atmfrc)
1516    end do
1517 #else
1518 
1519 ! Read by hand
1520    write(message,'(a,a,a,a)')'-Reading the file ',trim(filename),&
1521 & ' with Fortran'
1522 
1523    call wrtout(ab_out,message,'COLL')
1524    call wrtout(std_out,message,'COLL')
1525 
1526    if (open_file(filename,message,unit=funit,form="formatted",&
1527 &               status="old",action="read") /= 0) then
1528      MSG_ERROR(message)
1529    end if
1530 
1531 !Start a reading loop in fortran
1532    rewind(unit=funit)
1533    found=.false.
1534 
1535    iatom  = 1
1536    iamu   = 1
1537    itypat = 1
1538    irpt   = 1
1539    irpt1  = 0
1540    irpt2  = 0
1541    iph1l  = 1
1542    amu    = zero
1543    short_range  = .false.
1544    total_range  = .false.
1545 
1546    do while ((ios==0).and.(.not.found))
1547      read(funit,'(a)',iostat=ios) readline
1548      if(ios == 0)then
1549        call rmtabfromline(readline)
1550        line=adjustl(readline)
1551        if (.not.has_straincoupling) then
1552 
1553          if ((line(1:7)=='<energy')) then
1554            call rdfromline_value('energy',line,strg)
1555            if (strg/="") then
1556              strg1=trim(strg)
1557              read(strg1,*) energy
1558            else
1559              read(funit,'(a)',iostat=ios) readline
1560              call rmtabfromline(readline)
1561              line=adjustl(readline)
1562              call rdfromline_value('energy',line,strg)
1563              if (strg/="") then
1564                strg1=trim(strg)
1565              else
1566                strg1=trim(line)
1567              end if
1568              read(strg1,*) energy
1569            end if
1570            cycle
1571          end if
1572 
1573          if ((line(1:10)=='<unit_cell')) then
1574            call rdfromline_value('unit_cell',line,strg)
1575            if (strg/="") then
1576              strg1=trim(strg)
1577              read(strg1,*) (rprimd(1,mu),mu=1,3)
1578              read(funit,*) (rprimd(2,mu),mu=1,3)
1579            else
1580              do nu=1,2
1581                read(funit,*) (rprimd(nu,mu),mu=1,3)
1582              end do
1583            end if
1584            read(funit,'(a)',iostat=ios) readline
1585            call rmtabfromline(readline)
1586            line=adjustl(readline)
1587            call rdfromline_value('unit_cell',line,strg)
1588            if (strg/="") then
1589              strg1=trim(strg)
1590            else
1591              strg1=trim(line)
1592            end if
1593            read(strg1,*) (rprimd(3,mu),mu=1,3)
1594            cycle
1595          end if
1596 
1597          if ((line(1:12)=='<epsilon_inf')) then
1598            call rdfromline_value('epsilon_inf',line,strg)
1599            if (strg/="") then
1600              strg1=trim(strg)
1601              read(strg1,*) (epsilon_inf(mu,1),mu=1,3)
1602              read(funit,*) (epsilon_inf(mu,2),mu=1,3)
1603            else
1604              do nu=1,2
1605                read(funit,*) (epsilon_inf(mu,nu),mu=1,3)
1606              end do
1607            end if
1608            read(funit,'(a)',iostat=ios) readline
1609            call rmtabfromline(readline)
1610            line=adjustl(readline)
1611            call rdfromline_value('epsilon_inf',line,strg)
1612            if (strg/="") then
1613              strg1=trim(strg)
1614            else
1615              strg1=trim(line)
1616            end if
1617            read(strg1,*) (epsilon_inf(mu,3),mu=1,3)
1618            cycle
1619          end  if
1620 
1621          if ((line(1:8)=='<elastic')) then
1622            call rdfromline_value('elastic',line,strg)
1623            if (strg/="") then
1624              strg1=trim(strg)
1625              read(strg1,*) (elastic_constants(mu,1),mu=1,6)
1626              do nu=2,5
1627                read(funit,*) (elastic_constants(mu,nu),mu=1,6)
1628              end do
1629            else
1630              do nu=1,5
1631                read(funit,*) (elastic_constants(mu,nu),mu=1,6)
1632              end do
1633            end if
1634            read(funit,'(a)',iostat=ios) readline
1635            call rmtabfromline(readline)
1636            line=adjustl(readline)
1637            call rdfromline_value('elastic',line,strg)
1638            if (strg/="") then
1639              strg1=trim(strg)
1640            else
1641              strg1=trim(line)
1642            end if
1643            read(strg1,*) (elastic_constants(mu,6),mu=1,6)
1644            cycle
1645          end if
1646 
1647          if ((line(1:5)=='<atom')) then
1648            call rdfromline("mass",line,strg)
1649            strg1=trim(strg)
1650            read(strg1,*) amu
1651            if (.not.any(abs(all_amu-amu)<tol16)) then
1652              all_amu(iamu) = amu
1653              typat(iatom) = int(amu)
1654              !convert atomic mass unit to znucl
1655              do ii=1,103
1656                call atomdata_from_znucl(atom,real(ii,dp))
1657                if (abs((real(atom%amu,sp)-real(amu,sp))&
1658 &                 /real(amu,sp)*100)<0.1) then
1659                  znucl(iamu) = atom%znucl
1660                  exit
1661                end if
1662              end do
1663              iamu = iamu +1
1664            end if
1665            do itypat=1,ntypat
1666              if(abs(amu-all_amu(itypat))<tol16) then
1667                typat(iatom) = itypat
1668              end if
1669            end do
1670            cycle
1671          end if
1672 
1673          if ((line(1:9)=='<position')) then
1674            call rdfromline_value('position',line,strg)
1675            if (strg/="") then
1676              strg1=trim(strg)
1677              read(strg1,*)(xcart(mu,iatom),mu=1,3)
1678            else
1679              read(funit,'(a)',iostat=ios) readline
1680              call rmtabfromline(readline)
1681              line=adjustl(readline)
1682              call rdfromline_value('position',line,strg)
1683              if (strg/="") then
1684                strg1=trim(strg)
1685              else
1686                strg1=trim(line)
1687              end if
1688              read(strg1,*)(xcart(mu,iatom),mu=1,3)
1689            end if
1690            cycle
1691          end if
1692 
1693          if ((line(1:11)=='<borncharge')) then
1694            call rdfromline_value('borncharge',line,strg)
1695            if (strg/="") then
1696              strg1=trim(strg)
1697              read(strg1,*) (zeff(mu,1,iatom),mu=1,3)
1698              read(funit,*) (zeff(mu,2,iatom),mu=1,3)
1699            else
1700              do nu=1,2
1701                read(funit,*) (zeff(mu,nu,iatom),mu=1,3)
1702              end do
1703            end if
1704            read(funit,'(a)',iostat=ios) readline
1705            line=adjustl(readline)
1706            call rdfromline_value('borncharge',line,strg)
1707            if (strg/="") then
1708              strg1=trim(strg)
1709            else
1710              strg1=trim(line)
1711            end if
1712              read(strg1,*) (zeff(mu,3,iatom),mu=1,3)
1713            cycle
1714          end if
1715 
1716          if ((line(1:7)==char(60)//char(47)//'atom'//char(62))) then
1717            iatom=iatom+1
1718            cycle
1719          end if
1720 
1721          if ((line(1:12)=='<local_force')) then
1722            found2 = .False.
1723            irpt1 = irpt1 + 1
1724            do while (.not.found2)
1725              read(funit,'(a)',iostat=ios) readline
1726              call rmtabfromline(readline)
1727              line=adjustl(readline)
1728              if ((line(1:5)=='<data')) then
1729                call rdfromline_value('data',line,strg)
1730                if (strg/="") then
1731                  ABI_ALLOCATE(work2,(3*natom,3*natom))
1732                  strg1=trim(strg)
1733                  read(strg1,*) (work2(1,nu),nu=1,3*natom)
1734                  do mu=2,3*natom-1
1735                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
1736                  end do
1737                  read(funit,'(a)',iostat=ios) readline
1738                  call rmtabfromline(readline)
1739                  line=adjustl(readline)
1740                  call rdfromline_value('data',line,strg)
1741                  if (strg/="") then
1742                    strg1=trim(strg)
1743                  else
1744                    strg1=trim(line)
1745                  end if
1746                  read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1747                  local_atmfrc(:,:,:,:,irpt1) = reshape(work2,(/3,natom,3,natom/))
1748                  ABI_DEALLOCATE(work2)
1749                else
1750                  ABI_ALLOCATE(work2,(3*natom,3*natom))
1751                  do mu=1,3*natom
1752                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
1753                  end do
1754                  local_atmfrc(:,:,:,:,irpt1) =  reshape(work2,(/3,natom,3,natom/))
1755                  ABI_DEALLOCATE(work2)
1756                end if
1757              end if
1758              if ((line(1:5)=='<cell')) then
1759                call rdfromline_value('cell',line,strg)
1760                if (strg/="") then
1761                  strg1=trim(strg)
1762                  read(strg1,*)(cell_local(mu,irpt1),mu=1,3)
1763                else
1764                  read(funit,*)(cell_local(mu,irpt1),mu=1,3)
1765                end if
1766                found2 = .TRUE.
1767                cycle
1768              end if
1769            end do
1770          end if
1771 
1772          if ((line(1:12)=='<total_force')) then
1773            irpt2 = irpt2 + 1
1774            found2 = .False.
1775            do while (.not.found2)
1776              read(funit,'(a)',iostat=ios) readline
1777              call rmtabfromline(readline)
1778              line=adjustl(readline)
1779              if ((line(1:5)=='<data')) then
1780                call rdfromline_value('data',line,strg)
1781                if (strg/="") then
1782                  ABI_ALLOCATE(work2,(3*natom,3*natom))
1783                  strg1=trim(strg)
1784                  read(strg1,*) (work2(1,nu),nu=1,3*natom)
1785                  do mu=2,3*natom-1
1786                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
1787                  end do
1788                  read(funit,'(a)',iostat=ios) readline
1789                  call rmtabfromline(readline)
1790                  line=adjustl(readline)
1791                  call rdfromline_value('data',line,strg)
1792                  if (strg/="") then
1793                    strg1=trim(strg)
1794                  else
1795                    strg1=trim(line)
1796                  end if
1797                  read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
1798                  total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/))
1799                  ABI_DEALLOCATE(work2)
1800                else
1801                  ABI_ALLOCATE(work2,(3*natom,3*natom))
1802                  do mu=1,3*natom
1803                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
1804                  end do
1805                  total_atmfrc(:,:,:,:,irpt2) = reshape(work2,(/3,natom,3,natom/))
1806                  ABI_DEALLOCATE(work2)
1807                end if
1808              end if
1809              if ((line(1:5)=='<cell')) then
1810                call rdfromline_value('cell',line,strg)
1811                if (strg/="") then
1812                  strg1=trim(strg)
1813                  read(strg1,*)(cell_total(mu,irpt2),mu=1,3)
1814                else
1815                  read(funit,*)(cell_total(mu,irpt2),mu=1,3)
1816                end if
1817                found2 = .TRUE.
1818                cycle
1819              end if
1820            end do
1821          end if
1822 
1823          if ((line(1:7)=='<qpoint')) then
1824            call rdfromline_value('qpoint',line,strg)
1825            if (strg/="") then
1826              strg1=trim(strg)
1827              read(strg1,*)(qph1l(mu,iph1l),mu=1,3)
1828            else
1829              read(funit,*) (qph1l(mu,iph1l),mu=1,3)
1830            end if
1831          end if
1832 
1833          if ((line(1:12)=='<frequencies')) then
1834            call rdfromline_value('frequencies',line,strg)
1835            if (strg/="") then
1836              strg1=trim(strg)
1837              read(strg1,*)(phfrq(mu,iph1l),mu=1,3*natom)
1838            else
1839              do nu=1,natom
1840                read(funit,*) (phfrq(((nu-1)*3)+mu,iph1l),mu=1,3)
1841              end do
1842            end if
1843          end if
1844 
1845          if ((line(1:17)=='<dynamical_matrix')) then
1846            call rdfromline_value('dynamical_matrix',line,strg)
1847            if (strg/="") then
1848              ABI_ALLOCATE(work2,(3*natom,3*natom))
1849              strg1=trim(strg)
1850              read(strg1,*) (work2(nu,1),nu=1,3*natom)
1851              do mu=2,3*natom-1
1852                read(funit,*)(work2(nu,mu),nu=1,3*natom)
1853              end do
1854              read(funit,'(a)',iostat=ios) readline
1855              call rmtabfromline(readline)
1856              line=adjustl(readline)
1857              call rdfromline_value('dynamical_matrix',line,strg)
1858              if (strg/="") then
1859                strg1=trim(strg)
1860              else
1861                strg1=trim(line)
1862              end if
1863              read(strg1,*) (work2(nu,3*natom),nu=1,3*natom)
1864              dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/))
1865              ABI_DEALLOCATE(work2)
1866            else
1867              ABI_ALLOCATE(work2,(3*natom,3*natom))
1868              do mu=1,3*natom
1869                read(funit,*)(work2(nu,mu),nu=1,3*natom)
1870              end do
1871              dynmat(1,:,:,:,:,iph1l) = reshape(work2,(/3,natom,3,natom/))
1872              ABI_DEALLOCATE(work2)
1873            end if
1874          end if
1875 
1876          if ((line(1:8)==char(60)//char(47)//'phonon')) then
1877            iph1l = iph1l +1
1878          end if
1879 
1880          if ((line(1:16)=='<strain_coupling')) then
1881            read(funit,'(a)',iostat=ios) readline
1882            call rdfromline("voigt",line,strg)
1883            strg1=trim(strg)
1884            read(strg1,*) voigt
1885            voigt = voigt + 1 ! 0 to 5 in the xml
1886            has_straincoupling = .true.
1887            irpt = 1
1888          end if
1889 
1890        else
1891 !        Now treat the strain phonon coupling part
1892          if ((line(1:16)=='<strain_coupling')) then
1893            read(funit,'(a)',iostat=ios) readline
1894            call rdfromline("voigt",line,strg)
1895            strg1=trim(strg)
1896            read(strg1,*) voigt
1897            voigt = voigt + 1 ! 0 to 5 in the xml
1898            irpt = 1
1899            cycle
1900          end if
1901 
1902          if(voigt>6)then
1903            write(message, '(4a)' )ch10,&
1904 &               ' WARNING: the number of strain phonon coupling is superior to 6 in ',filename,ch10
1905            call wrtout(std_out,message,'COLL')
1906            exit
1907          end if
1908 
1909          if ((line(1:22)=='<correction_force unit')) then
1910            call rdfromline_value('correction_force',line,strg)
1911            if (strg/="") then
1912              ABI_ALLOCATE(work2,(3,natom))
1913              strg1=trim(strg)
1914              read(strg1,*) (work2(nu,1),nu=1,3)
1915              do mu=2,natom-1
1916                read(funit,*)(work2(nu,mu),nu=1,3)
1917              end do
1918              read(funit,'(a)',iostat=ios) readline
1919              call rmtabfromline(readline)
1920              line=adjustl(readline)
1921              call rdfromline_value('correction_force',line,strg)
1922              if (strg/="") then
1923                strg1=trim(strg)
1924                read(strg1,*) (work2(nu,natom),nu=1,3)
1925              else
1926                strg1=trim(line)
1927                read(strg1,*) (work2(nu,natom),nu=1,3)
1928              end if
1929              strain_coupling(voigt,:,:) = work2(:,:)
1930              ABI_DEALLOCATE(work2)
1931            else
1932              ABI_ALLOCATE(work2,(3,natom))
1933              do mu=1,natom
1934                read(funit,*)(work2(nu,mu),nu=1,3)
1935              end do
1936              strain_coupling(voigt,:,:) = work2(:,:)
1937              ABI_DEALLOCATE(work2)
1938            end if
1939          end if
1940 
1941          if ((line(1:11)=='<elastic3rd')) then
1942            call rdfromline_value('elastic3rd',line,strg)
1943            if (strg/="") then
1944              strg1=trim(strg)
1945              read(strg1,*) (elastic3rd(voigt,mu,1),mu=1,6)
1946              do nu=2,5
1947                read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6)
1948              end do
1949            else
1950              do nu=1,5
1951                read(funit,*) (elastic3rd(voigt,mu,nu),mu=1,6)
1952              end do
1953            end if
1954            read(funit,'(a)',iostat=ios) readline
1955            call rmtabfromline(readline)
1956            line=adjustl(readline)
1957            call rdfromline_value('elastic3rd',line,strg)
1958            if (strg/="") then
1959              strg1=trim(strg)
1960              read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6)
1961            else
1962              strg1=trim(line)
1963              read(strg1,*) (elastic3rd(voigt,mu,6),mu=1,6)
1964            end if
1965            has_anharmonics = .true.
1966            cycle
1967          end if
1968 
1969          if ((line(1:29)=='<correction_strain_force unit')) then
1970            call rdfromline_value('correction_strain_force',line,strg)
1971            if (strg/="") then
1972              ABI_ALLOCATE(work2,(3*6,natom))
1973              strg1=trim(strg)
1974              read(strg1,*) (work2(nu,1),nu=1,3*6)
1975              do mu=2,natom-1
1976                read(funit,*)(work2(nu,mu),nu=1,3*6)
1977              end do
1978              read(funit,'(a)',iostat=ios) readline
1979              call rmtabfromline(readline)
1980              line=adjustl(readline)
1981              call rdfromline_value('correction_strain_force',line,strg)
1982              if (strg/="") then
1983                strg1=trim(strg)
1984                read(strg1,*) (work2(nu,natom),nu=1,3*6)
1985              else
1986                strg1=trim(line)
1987                read(strg1,*) (work2(nu,natom),nu=1,3*6)
1988              end if
1989              elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/))
1990              ABI_DEALLOCATE(work2)
1991            else
1992              ABI_ALLOCATE(work2,(3*6,natom))
1993              do mu=1,natom
1994                read(funit,*)(work2(nu,mu),nu=1,3*6)
1995              end do
1996              elastic_displacement(voigt,:,:,:) = reshape(work2(:,:),(/6,3,natom/))
1997              ABI_DEALLOCATE(work2)
1998            end if
1999          end if
2000 
2001          if ((line(1:26)=='<correction_force_constant')) then
2002            found2=.false.
2003            do while (.not.found2)
2004              read(funit,'(a)',iostat=ios) readline
2005              call rmtabfromline(readline)
2006              line=adjustl(readline)
2007              if ((line(1:5)=='<data')) then
2008                call rdfromline_value('data',line,strg)
2009                if (strg/="") then
2010                  ABI_ALLOCATE(work2,(3*natom,3*natom))
2011                  strg1=trim(strg)
2012                  read(strg1,*) (work2(1,nu),nu=1,3*natom)
2013                  do mu=2,3*natom-1
2014                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
2015                  end do
2016                  read(funit,'(a)',iostat=ios) readline
2017                  call rmtabfromline(readline)
2018                  line=adjustl(readline)
2019                  call rdfromline_value('data',line,strg)
2020                  if (strg/="") then
2021                    strg1=trim(strg)
2022                    read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
2023                  else
2024                    strg1=trim(line)
2025                    read(strg1,*) (work2(3*natom,nu),nu=1,3*natom)
2026                  end if
2027                  phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) = &
2028 &                           reshape(work2,(/3,natom,3,natom/))
2029                  ABI_DEALLOCATE(work2)
2030                else
2031                  ABI_ALLOCATE(work2,(3*natom,3*natom))
2032                  do mu=1,3*natom
2033                    read(funit,*)(work2(mu,nu),nu=1,3*natom)
2034                  end do
2035                  phonon_strain(voigt)%atmfrc(:,:,:,:,irpt) =&
2036 &              reshape(work2,(/3,natom,3,natom/))
2037                  ABI_DEALLOCATE(work2)
2038                end if
2039                has_anharmonics = .true.
2040              end if
2041              if ((line(1:5)=='<cell')) then
2042                call rdfromline_value('cell',line,strg)
2043                if (strg/="") then
2044                  strg1=trim(strg)
2045                  read(strg1,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3)
2046                else
2047                  read(funit,*)(phonon_strain(voigt)%cell(mu,irpt),mu=1,3)
2048                end if
2049                irpt = irpt + 1
2050                found2=.true.
2051                cycle
2052              end if
2053            end do
2054          end if
2055 
2056          if ((line(1:17)==char(60)//char(47)//'strain_coupling')) then
2057 !          set nrpt for the previous value of strain
2058            phonon_strain(voigt)%nrpt = irpt - 1
2059 !         restart the calculation of nrpt
2060          end if
2061        end if
2062      end if
2063    end do
2064 
2065 
2066 ! Reorder the ATMFRC
2067 ! Case 1: only local in the xml
2068    if (irpt1>0 .and. irpt2==0) then
2069      ifcs%cell(:,:) = int(cell_local(:,:))
2070      ifcs%atmfrc(:,:,:,:,:)  = local_atmfrc(:,:,:,:,:)
2071      ifcs%short_atmfrc(:,:,:,:,:) = local_atmfrc(:,:,:,:,:)
2072      ifcs%ewald_atmfrc(:,:,:,:,:) = zero
2073 
2074 ! Case 2: only total in the xml
2075    else if(irpt1==0 .and. irpt2>0)then
2076      ifcs%cell(:,:) = int(cell_total(:,:))
2077      ifcs%atmfrc(:,:,:,:,:)  = total_atmfrc(:,:,:,:,:)
2078      ifcs%short_atmfrc(:,:,:,:,:) = zero
2079      ifcs%ewald_atmfrc(:,:,:,:,:) = total_atmfrc(:,:,:,:,:)
2080 
2081 ! Case 3: local + total in the xml
2082    else if (irpt1>0 .and. irpt2>0)then
2083      if(irpt1 <= irpt2)then
2084        irpt3 = 0
2085        do ii=1,irpt2
2086          ifcs%cell(:,ii) = int(cell_total(:,ii))
2087          ifcs%atmfrc(:,:,:,:,ii)  = total_atmfrc(:,:,:,:,ii)
2088          do jj=1,irpt1
2089            if (all(abs(int(cell_local(:,jj))-ifcs%cell(:,ii))<tol16)) then
2090              ifcs%short_atmfrc(:,:,:,:,ii) = local_atmfrc(:,:,:,:,jj)
2091              irpt3 = irpt3 + 1
2092            end if
2093          end do
2094        end do
2095        if(irpt3 /= irpt1)then
2096          write(message, '(4a)' )ch10,&
2097 &         ' There is several similar short IFC in ',filename,ch10
2098          MSG_BUG(message)
2099        end if
2100      else
2101        write(message, '(2a,I5,3a,I5,5a)' )ch10,&
2102 &     ' The number of total IFC  (',irpt2,') is inferior to  ',ch10,&
2103 &     ' the number of short range IFC (',irpt1,') in ',filename,ch10,&
2104 &     ' This is not possible',ch10
2105        MSG_BUG(message)
2106      end if
2107    end if
2108 
2109 !  Do some checks
2110    if (any(typat==0)) then
2111      write(message, '(a,a,a)' )&
2112 &      ' Unable to read the type of atoms ',trim(filename),ch10
2113      MSG_ERROR(message)
2114    end if
2115 
2116    if (any(abs(znucl)<tol16)) then
2117      write(message, '(a,a,a)' )&
2118 &      ' Unable to read the atomic number ',trim(filename),ch10
2119      MSG_ERROR(message)
2120    end if
2121 
2122    if (any(abs(all_amu)<tol16)) then
2123      write(message, '(a,a,a)' )&
2124 &     ' Unable to read the atomic mass ',trim(filename),ch10
2125      MSG_ERROR(message)
2126    end if
2127 
2128    close(unit=funit)
2129 
2130 #endif
2131 
2132  end if !End if master
2133 
2134 !MPI BROADCAST
2135  call xmpi_bcast(energy,master, comm, ierr)
2136  call xmpi_bcast(all_amu,master, comm, ierr)
2137  call xmpi_bcast(dynmat,master, comm, ierr)
2138  call xmpi_bcast(elastic_constants,master, comm, ierr)
2139  call xmpi_bcast(epsilon_inf,master, comm, ierr)
2140  call xmpi_bcast(ifcs%nrpt,master, comm, ierr)
2141  call xmpi_bcast(ifcs%atmfrc,master, comm, ierr)
2142  call xmpi_bcast(ifcs%cell,master, comm, ierr)
2143  call xmpi_bcast(ifcs%ewald_atmfrc,master, comm, ierr)
2144  call xmpi_bcast(ifcs%short_atmfrc,master, comm, ierr)
2145  call xmpi_bcast(strain_coupling,master, comm, ierr)
2146  call xmpi_bcast(phfrq,master, comm, ierr)
2147  call xmpi_bcast(qph1l,master, comm, ierr)
2148  call xmpi_bcast(typat,master, comm, ierr)
2149  call xmpi_bcast(rprimd,master, comm, ierr)
2150  call xmpi_bcast(xcart,master, comm, ierr)
2151  call xmpi_bcast(zeff,master, comm, ierr)
2152  call xmpi_bcast(znucl,master, comm, ierr)
2153  do ii = 1,6
2154    call xmpi_bcast(phonon_strain(ii)%nrpt   ,master, comm, ierr)
2155    call xmpi_bcast(phonon_strain(ii)%atmfrc ,master, comm, ierr)
2156    call xmpi_bcast(phonon_strain(ii)%cell   ,master, comm, ierr)
2157  end do
2158  call xmpi_bcast(elastic3rd   ,master, comm, ierr)
2159  call xmpi_bcast(elastic_displacement ,master, comm, ierr)
2160  call xmpi_bcast(has_anharmonics ,master, comm, ierr)
2161 
2162 !Fill somes others variables
2163  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2164  call xcart2xred(natom,rprimd,xcart,xred)
2165 
2166 !Re-generate symmetry operations from the lattice and atomic coordinates
2167  tolsym=tol8
2168  msym = 192
2169  ABI_ALLOCATE(spinat,(3,natom))
2170  ABI_ALLOCATE(ptsymrel,(3,3,msym))
2171  ABI_ALLOCATE(symafm,(msym))
2172  ABI_ALLOCATE(symrel,(3,3,msym))
2173  ABI_ALLOCATE(tnons,(3,msym))
2174  use_inversion=1
2175  spinat = 0;
2176  symrel = 0;
2177  symafm = 0;
2178  tnons = 0 ;
2179  space_group = 0;
2180  call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
2181  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2182  call symfind(0,(/zero,zero,zero/),gprimd,0,msym,natom,0,nptsym,nsym,&
2183 &  0,0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred)
2184 
2185 !Initialisation of crystal
2186  npsp = ntypat; timrev = 1
2187  ABI_ALLOCATE(title, (ntypat))
2188  do ii=1,ntypat
2189    write(title(ii),'(a,i0)')"No title for typat ",ii
2190  end do
2191 
2192 !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
2193  call crystal_init(all_amu,Crystal,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,&
2194 &  zion,znucl,timrev,.FALSE.,.FALSE.,title,&
2195 &  symrel=symrel(:,:,1:nsym),tnons=tnons(:,1:nsym),symafm=symafm(1:nsym))
2196 
2197 !amu is not fill in crystal_init...
2198  Crystal%amu(:) = all_amu(:)
2199 
2200  ABI_DEALLOCATE(symrel)
2201  ABI_DEALLOCATE(symafm)
2202  ABI_DEALLOCATE(tnons)
2203  ABI_DEALLOCATE(spinat)
2204  ABI_DEALLOCATE(ptsymrel)
2205 
2206 !if strcpling is set to 0 by the user, need to set the flag to false for
2207 !the initialisation of the effective potential
2208  if (present(strcpling))then
2209    if(strcpling == 0 )then
2210      has_anharmonics = .FALSE.
2211    end if
2212  end if
2213 
2214 !Initialisation of eff_pot
2215  call effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nph1l,comm,&
2216 &                              dynmat=dynmat,elastic_constants=elastic_constants,&
2217 &                              elastic3rd=elastic3rd,elastic_displacement=elastic_displacement,&
2218 &                              epsilon_inf=epsilon_inf,strain_coupling=strain_coupling,&
2219 &                              phonon_strain=phonon_strain,phfrq=phfrq,qpoints=qph1l,&
2220 &                              has_anharmonicsTerms=has_anharmonics,zeff=zeff)
2221 
2222 !DEALLOCATION OF ARRAYS
2223  ABI_DEALLOCATE(all_amu)
2224  ABI_DEALLOCATE(cell_local)
2225  ABI_DEALLOCATE(cell_total)
2226  ABI_DEALLOCATE(total_atmfrc)
2227  ABI_DEALLOCATE(local_atmfrc)
2228  ABI_DEALLOCATE(ifcs%atmfrc)
2229  ABI_DEALLOCATE(ifcs%cell)
2230  ABI_DEALLOCATE(ifcs%short_atmfrc)
2231  ABI_DEALLOCATE(ifcs%ewald_atmfrc)
2232  ABI_DEALLOCATE(dynmat)
2233  ABI_DEALLOCATE(strain_coupling)
2234  ABI_DEALLOCATE(phfrq)
2235  ABI_DEALLOCATE(qph1l)
2236  ABI_DEALLOCATE(title)
2237  ABI_DEALLOCATE(typat)
2238  ABI_DEALLOCATE(xcart)
2239  ABI_DEALLOCATE(xred)
2240  ABI_DEALLOCATE(zeff)
2241  ABI_DEALLOCATE(zion)
2242  ABI_DEALLOCATE(znucl)
2243  do ii = 1,6
2244    phonon_strain(ii)%nrpt   = nrpt
2245    phonon_strain(ii)%atmfrc = zero
2246    phonon_strain(ii)%cell   = 0
2247    ABI_DEALLOCATE(phonon_strain(ii)%atmfrc)
2248    ABI_DEALLOCATE(phonon_strain(ii)%cell)
2249  end do
2250  ABI_FREE(phonon_strain)
2251  ABI_DEALLOCATE(elastic_displacement)
2252 
2253 !DEALLOCATION OF TYPES
2254  call ifc_free(ifcs)
2255  call crystal_free(crystal)
2256 
2257 end subroutine system_xml2effpot

m_effpot_xml/char_c2f [ Functions ]

[ Top ] [ Functions ]

NAME

  char_c_to_f

FUNCTION

 Helper function to convert a C string to a Fortran string
 Based on a routine by Joseph M. Krahn

INPUTS

  c_string=C string

OUTPUT

  f_string=Fortran string

PARENTS

      m_libpaw_libxc

CHILDREN

SOURCE

4111 subroutine char_c2f(c_string,f_string)
4112 
4113  use iso_c_binding, only : C_CHAR,C_NULL_CHAR
4114 !Arguments ------------------------------------
4115 
4116 !This section has been created automatically by the script Abilint (TD).
4117 !Do not modify the following lines by hand.
4118 #undef ABI_FUNC
4119 #define ABI_FUNC 'char_c2f'
4120 !End of the abilint section
4121 
4122  character(kind=C_CHAR,len=1),intent(in) :: c_string(*)
4123  character(len=*),intent(out) :: f_string
4124 !Local variables -------------------------------
4125  integer :: ii
4126 !! *************************************************************************
4127  ii=1
4128  do while(c_string(ii)/=C_NULL_CHAR.and.ii<=len(f_string))
4129    f_string(ii:ii)=c_string(ii) ; ii=ii+1
4130  end do
4131  if (ii<len(f_string)) f_string(ii:)=' '
4132 end subroutine char_c2f

m_effpot_xml/char_f2c [ Functions ]

[ Top ] [ Functions ]

NAME

  char_f_to_c

FUNCTION

 Helper function to convert a Fortran string to a C string
 Based on a routine by Joseph M. Krahn

INPUTS

  f_string=Fortran string

OUTPUT

  c_string=C string

SOURCE

4062 #if defined HAVE_LIBXML
4063 
4064 function char_f2c(f_string) result(c_string)
4065 
4066  use iso_c_binding, only : C_CHAR,C_NULL_CHAR
4067 !Arguments ------------------------------------
4068 
4069 !This section has been created automatically by the script Abilint (TD).
4070 !Do not modify the following lines by hand.
4071 #undef ABI_FUNC
4072 #define ABI_FUNC 'char_f2c'
4073 !End of the abilint section
4074 
4075  character(len=*),intent(in) :: f_string
4076  character(kind=C_CHAR,len=1) :: c_string(len_trim(f_string)+1)
4077 !Local variables -------------------------------
4078  integer :: ii,strlen
4079 !! *************************************************************************
4080  strlen=len_trim(f_string)
4081  forall(ii=1:strlen)
4082    c_string(ii)=f_string(ii:ii)
4083  end forall
4084  c_string(strlen+1)=C_NULL_CHAR
4085 end function char_f2c