TABLE OF CONTENTS


ABINIT/m_fit_polynomial_coeff [ Modules ]

[ Top ] [ Modules ]

NAME

 m_fit_polynomial_coeff

FUNCTION

COPYRIGHT

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

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 module m_fit_polynomial_coeff
26 
27  use defs_basis
28  use m_errors
29  use m_abicore
30  use m_polynomial_coeff
31  use m_atomdata
32  use m_xmpi
33  use m_supercell
34 
35  use m_special_funcs,only : factorial
36  use m_geometry,       only : xred2xcart
37  use m_crystal,only : symbols_crystal
38  use m_strain,only : strain_type,strain_get
39  use m_effective_potential,only : effective_potential_type, effective_potential_evaluate
40  use m_effective_potential,only : effective_potential_freeCoeffs,effective_potential_setCoeffs
41  use m_effective_potential,only : effective_potential_getDisp
42  use m_effective_potential_file, only : effective_potential_file_mapHistToRef
43  use m_io_tools,   only : open_file,get_unit
44  use m_abihist, only : abihist,abihist_free,abihist_init,abihist_copy,write_md_hist,var2hist
45  use m_random_zbq
46  use m_fit_data
47 
48  implicit none
49 
50  public :: fit_polynomial_coeff_computeGF
51  public :: fit_polynomial_coeff_computeMSD
52  public :: fit_polynomial_coeff_fit
53  public :: fit_polynomial_coeff_getFS
54  public :: fit_polynomial_coeff_getPositive
55  public :: fit_polynomial_coeff_getCoeffBound
56  public :: fit_polynomial_coeff_solve
57  public :: fit_polynomial_printSystemFiles
58  public :: genereList

m_fit_polynomial_coeff/fit_polynomial_coeff_computeGF [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_computeGF

FUNCTION

 Compute the values of the goal function (Mean squared error) for
   gf_value(1) = forces (Ha/Bohr)**2
   gf_value(2) = stresses (Ha/Bohr)**2
   gf_value(3) = stresses+forces (Ha/Bohr)**2
   gf_value(4) = energy (Ha)

INPUTS

 coefficients(ncoeff)          = type(polynomial_coeff_type)
 energy_coeffs(ncoeff,ntime)   = value of the energy for each  coefficient (Ha)
 energy_diff(ntime) = Difference of energ ybetween DFT calculation and fixed part
                             of the model (more often harmonic part)
                             fixed part of the model (more often harmonic part)
 fcart_coeffs(ncoeff,3,natom,ntime) = value of the forces for each coefficient
                                      (-1 factor is taking into acount) (Ha/Bohr)
 fcart_diff(3,natom,ntime) = Difference of cartesian forces between DFT calculation and
                             fixed part of the model (more often harmonic part)
 list_coeffs(ncoeff_fit) = List with the indexes of the coefficients used for this model
 natom = Number of atoms
 ncoeff_fit = Number of coefficients fitted
 ncoeff_max = Maximum number of coeff in the list
 ntime = Number of time in the history
 strten_coeffs(ncoeff,3,natom,ntime)= value of the stresses for each coefficient
                                      (1/ucvol factor is taking into acount) (Ha/Bohr^3)
 strten_diff(6,natom) = Difference of stress tensor between DFT calculation and
                        fixed part of the model (more often harmonic part)
 sqomega =  Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]

OUTPUT

 gf_value(4) = Goal function

PARENTS

      m_fit_polynomial_coeff

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

1838 subroutine fit_polynomial_coeff_computeGF(coefficients,energy_coeffs,energy_diff,&
1839 &                                         fcart_coeffs,fcart_diff,gf_value,list_coeffs,&
1840 &                                         natom,ncoeff_fit,ncoeff_max,ntime,strten_coeffs,&
1841 &                                         strten_diff,sqomega)
1842 
1843 
1844 !This section has been created automatically by the script Abilint (TD).
1845 !Do not modify the following lines by hand.
1846 #undef ABI_FUNC
1847 #define ABI_FUNC 'fit_polynomial_coeff_computeGF'
1848 !End of the abilint section
1849 
1850  implicit none
1851 
1852 !Arguments ------------------------------------
1853 !scalars
1854  integer,intent(in)  :: natom,ncoeff_fit,ncoeff_max,ntime
1855 !arrays
1856  integer,intent(in)  :: list_coeffs(ncoeff_fit)
1857  real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime)
1858  real(dp),intent(in) :: energy_diff(ntime)
1859  real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime)
1860  real(dp),intent(in) :: fcart_diff(3,natom,ntime)
1861  real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max)
1862  real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime)
1863  real(dp),intent(in) :: coefficients(ncoeff_fit)
1864  real(dp),intent(out) :: gf_value(4)
1865 !Local variables-------------------------------
1866 !scalar
1867  integer :: ia,icoeff,icoeff_tmp,itime,mu
1868  real(dp):: etmp,emu,fmu,ftmp,smu,stmp
1869  real(dp) :: ffact,sfact,efact
1870 !arrays
1871 ! *************************************************************************
1872 
1873 !1-Compute the value of the goal function
1874 ! see equation 9 of PRB 95 094115(2017) [[cite:Escorihuela-Sayalero2017]]
1875  gf_value = zero
1876  etmp     = zero
1877  ftmp     = zero
1878  stmp     = zero
1879 
1880 !Compute factors
1881  ffact = one/(3*natom*ntime)
1882  sfact = one/(6*ntime)
1883  efact = one/(ntime)
1884  
1885 ! loop over the configuration
1886  do itime=1,ntime
1887 ! Fill energy
1888    emu = zero
1889    do icoeff=1,ncoeff_fit
1890      icoeff_tmp = list_coeffs(icoeff)
1891      emu = emu + coefficients(icoeff)*energy_coeffs(icoeff_tmp,itime)
1892    end do
1893 !   uncomment the next line to be consistent with the definition of the goal function   
1894 !   etmp = etmp + (energy_diff(itime)-emu)**2
1895    etmp = etmp + abs(energy_diff(itime)-emu)
1896 !  Fill forces
1897    do ia=1,natom
1898      do mu=1,3
1899        fmu  = zero
1900        do icoeff=1,ncoeff_fit
1901          icoeff_tmp = list_coeffs(icoeff)
1902          fmu =  fmu + coefficients(icoeff)*fcart_coeffs(mu,ia,icoeff_tmp,itime)
1903        end do
1904        ftmp = ftmp + (fcart_diff(mu,ia,itime)-fmu)**2
1905      end do !End loop dir
1906    end do !End loop natom
1907    do mu=1,6
1908      smu = zero
1909      do icoeff=1,ncoeff_fit
1910        icoeff_tmp = list_coeffs(icoeff)
1911        smu = smu + coefficients(icoeff)*strten_coeffs(mu,itime,icoeff_tmp)
1912      end do
1913      stmp = stmp + sqomega(itime)*(strten_diff(mu,itime)-smu)**2
1914    end do !End loop stress dir
1915  end do ! End loop time
1916 
1917  gf_value(1)   =  ffact*ftmp + sfact*stmp !+ efact*etmp !Stresses + Forces
1918  gf_value(2)   =  ffact*ftmp ! only Forces
1919  gf_value(3)   =  sfact*stmp ! only Stresses
1920  gf_value(4)   =  efact*etmp !abs(Energy)
1921 
1922 end subroutine fit_polynomial_coeff_computeGF

m_fit_polynomial_coeff/fit_polynomial_coeff_computeMSD [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_computeMSD

FUNCTION

 Compute the Mean square error of the energy, forces and stresses

INPUTS

 eff_pot<type(effective_potential)> = effective potential
 hist<type(abihist)> = The history of the MD
 natom = number of atom
 ntime = number of time in the hist
 sqomega =  Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]
 compute_anharmonic = TRUE if the anharmonic part of the effective potential
                           has to be taking into acount
 print_file = if True, a ASCII file with the difference in energy will be print

OUTPUT

 mse  =  Mean square error of the energy   (Hatree)
 msef =  Mean square error of the forces   (Hatree/Bohr)**2
 mses =  Mean square error of the stresses (Hatree/Bohr)**2

PARENTS

      m_fit_polynomial_coeff

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

2254 subroutine fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,sqomega,&
2255 &                                          compute_anharmonic,print_file)
2256 
2257 
2258 !This section has been created automatically by the script Abilint (TD).
2259 !Do not modify the following lines by hand.
2260 #undef ABI_FUNC
2261 #define ABI_FUNC 'fit_polynomial_coeff_computeMSD'
2262 !End of the abilint section
2263 
2264  implicit none
2265 
2266 !Arguments ------------------------------------
2267 !scalars
2268  integer, intent(in) :: natom,ntime
2269  real(dp),intent(out):: mse,msef,mses
2270  logical,optional,intent(in) :: compute_anharmonic,print_file
2271 !arrays
2272  real(dp) :: sqomega(ntime)
2273  type(effective_potential_type),intent(in) :: eff_pot
2274  type(abihist),intent(in) :: hist
2275 !Local variables-------------------------------
2276 !scalar
2277  integer :: ii,ia,mu,unit_energy,unit_stress
2278 ! integer :: ifirst
2279  real(dp):: energy,energy_harm
2280  logical :: need_anharmonic = .TRUE.,need_print=.FALSE.
2281  !arrays
2282  real(dp):: fcart(3,natom),fred(3,natom),strten(6),rprimd(3,3),xred(3,natom)
2283  character(len=500) :: msg
2284 ! type(abihist) :: hist_out
2285 ! character(len=200) :: filename_hist
2286 
2287 ! *************************************************************************
2288 
2289  !Do some checks
2290  if(ntime /= hist%mxhist)then
2291    write(msg,'(a)')'ntime is not correct'
2292    MSG_BUG(msg)
2293  end if
2294 
2295  if(natom /= size(hist%xred,2)) then
2296    write(msg,'(a)')'natom is not correct'
2297    MSG_BUG(msg)
2298  end if
2299 
2300  if(present(compute_anharmonic))then
2301    need_anharmonic = compute_anharmonic
2302  end if
2303 
2304  if(present(print_file))then
2305 !   call abihist_init(hist_out,natom,ntime,.false.,.false.)
2306    need_print=print_file
2307    unit_energy = get_unit()
2308    if (open_file('fit_diff_energy.dat',msg,unit=unit_energy,form="formatted",&
2309 &     status="unknown",action="write") /= 0) then
2310      MSG_ERROR(msg)
2311    end if
2312 
2313    unit_stress = get_unit()
2314    if (open_file('fit_diff_stress.dat',msg,unit=unit_stress,form="formatted",&
2315 &     status="unknown",action="write") /= 0) then
2316      MSG_ERROR(msg)
2317    end if
2318 
2319  end if
2320 
2321  mse  = zero
2322  msef = zero
2323  mses = zero
2324 
2325  do ii=1,ntime
2326    xred(:,:)   = hist%xred(:,:,ii)
2327    rprimd(:,:) = hist%rprimd(:,:,ii)
2328    call effective_potential_evaluate(eff_pot,energy_harm,fcart,fred,strten,natom,rprimd,&
2329 &                                    xred=xred,compute_anharmonic=.False.,verbose=.false.)
2330 
2331    call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd,&
2332 &                                    xred=xred,compute_anharmonic=need_anharmonic,verbose=.false.)
2333 
2334    if(need_print)then
2335      WRITE(unit_energy ,'(I10,5(F23.14))') ii,hist%etot(ii),energy_harm,energy,&
2336 &                                       abs(hist%etot(ii) - energy_harm),abs(hist%etot(ii) - energy)
2337      WRITE(unit_stress,'(I10,12(F23.14))') ii,hist%strten(:,ii),strten(:)
2338    end if
2339 
2340 !    ifirst=merge(0,1,(ii>1))
2341 !    filename_hist = trim("test.nc")
2342 !    hist_out%fcart(:,:,hist_out%ihist) = hist%fcart(:,:,ii)
2343 !    hist_out%strten(:,hist_out%ihist)  = hist%strten(:,ii)
2344 !    hist_out%etot(hist_out%ihist)      = hist%etot(ii)
2345 !    hist_out%entropy(hist_out%ihist)   = hist%entropy(ii)
2346 !    hist_out%time(hist_out%ihist)      = real(ii,kind=dp)
2347 !    call vel2hist(ab_mover%amass,hist,vel,vel_cell)
2348 !    call var2hist(hist%acell(:,ii),hist_out,natom,hist%rprimd(:,:,ii),hist%xred(:,:,ii),.false.)
2349 !    call write_md_hist(hist_out,filename_hist,ifirst,ii,natom,1,eff_pot%crystal%ntypat,&
2350 ! &                    eff_pot%supercell%typat,eff_pot%crystal%amu,eff_pot%crystal%znucl,&
2351 ! &                    real(100,dp),(/real(100,dp),real(100,dp)/))
2352 
2353    mse  = mse  + abs(hist%etot(ii) - energy)
2354    do ia=1,natom
2355      do mu=1,3
2356        msef = msef + (hist%fcart(mu,ia,ii)  - fcart(mu,ia))**2
2357      end do
2358    end do
2359    do mu=1,6
2360      mses = mses + sqomega(ii)*(hist%strten(mu,ii) - strten(mu))**2
2361    end do
2362  end do
2363 
2364  mse  = mse  /  ntime
2365  msef = msef / (3*natom*ntime)
2366  mses = mses / (6*ntime)
2367 
2368  if(need_print)then
2369    close(unit_energy)
2370    close(unit_stress)
2371  end if
2372 
2373 ! call abihist_free(hist_out)
2374 
2375 end subroutine fit_polynomial_coeff_computeMSD

m_fit_polynomial_coeff/fit_polynomial_coeff_fit [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_fit

FUNCTION

 Fit the list of coefficients included in eff_pot,
 if the coefficients are not set in eff_pot, this routine will genenerate
 a list of coefficients by taking into acount the symmetries of the system
 and the cutoff

INPUTS

 eff_pot<type(effective_potential)> = effective potential
 bancoeff(nbancoeff) = list of bannned coeffcients, these coefficients will NOT be
                       used during the fit process
 fixcoeff(nfixcoeff) = list of fixed coefficient, these coefficients will be
                       imposed during the fit process
 hist<type(abihist)> = The history of the MD (or snapshot of DFT)
 generateterm = term to activate the generation of the term set
 power_disps(2) = array with the minimal and maximal power_disp to be computed
 nbancoeff = number of banned coeffcients
 ncycle_in = number of maximum cycle (maximum coefficient to be fitted)
 nfixcoeff = Number of coefficients imposed during the fit process
 option = option of the fit process : 1 - selection of the coefficient one by one
                                      2 - selection of the coefficients with Monte Carlo(testversion)
 comm = MPI communicator
 cutoff_in = optional,cut off to apply to the range of interation if
           the coefficient are genereted in this routine
 max_power_strain = maximum order of the strain of the strain phonon coupling
 fit_initializeData = optional, logical !If true, we store all the informations for the fit,
                      it will reduce the computation time but increase a lot the memory...
 fit_tolMSDF = optional, tolerance in eV^2/A^2 on the Forces for the fit process
 fit_tolMSDS = optional, tolerance in eV^2/A^2 on the Stresses for the fit process
 fit_tolMSDE = optional, tolerance in meV^2/A^2 on the Energy for the fit process
 fit_tolMSDFS= optional, tolerance in eV^2/A^2 on the Forces+stresses for the fit process
 positive = optional, TRUE will return only positive coefficients
                      FALSE, default
 verbose  = optional, flag for the verbose mode
 anhstr = logical, optional : TRUE, the anharmonic strain are computed
                              FALSE, (default) the anharmonic strain are not computed
 only_odd_power = logical, optional : if TRUE generate only odd power
 only_even_power= logical, optional : if TRUE generate only even power

OUTPUT

 eff_pot<type(effective_potential)> = effective potential datatype with new fitted coefficients

PARENTS

      m_fit_polynomial_coeff,mover_effpot,multibinit

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

 118 subroutine fit_polynomial_coeff_fit(eff_pot,bancoeff,fixcoeff,hist,generateterm,power_disps,&
 119 &                                   nbancoeff,ncycle_in,nfixcoeff,option,comm,cutoff_in,&
 120 &                                   max_power_strain,initialize_data,&
 121 &                                   fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS,&
 122 &                                   positive,verbose,anharmstr,spcoupling,&
 123 &                                   only_odd_power,only_even_power)
 124 
 125 
 126 !This section has been created automatically by the script Abilint (TD).
 127 !Do not modify the following lines by hand.
 128 #undef ABI_FUNC
 129 #define ABI_FUNC 'fit_polynomial_coeff_fit'
 130 !End of the abilint section
 131 
 132  implicit none
 133 
 134 !Arguments ------------------------------------
 135 !scalars
 136  integer,intent(in) :: ncycle_in,nfixcoeff,comm
 137  integer,intent(in) :: generateterm,nbancoeff,option
 138 !arrays
 139  integer,intent(in) :: fixcoeff(nfixcoeff), bancoeff(nbancoeff)
 140  integer,intent(in) :: power_disps(2)
 141  type(effective_potential_type),target,intent(inout) :: eff_pot
 142  type(abihist),intent(inout) :: hist
 143  integer,optional,intent(in) :: max_power_strain
 144  real(dp),optional,intent(in) :: cutoff_in,fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS
 145  logical,optional,intent(in) :: verbose,positive,anharmstr,spcoupling
 146  logical,optional,intent(in) :: only_odd_power,only_even_power
 147  logical,optional,intent(in) :: initialize_data
 148 !Local variables-------------------------------
 149 !scalar
 150  integer :: ii,icoeff,my_icoeff,icycle,icycle_tmp,ierr,info,index_min,iproc,isweep,jcoeff
 151  integer :: master,max_power_strain_in,my_rank,my_ncoeff,ncoeff_model,ncoeff_tot,natom_sc,ncell,ncycle
 152  integer :: ncycle_tot,ncycle_max,nproc,ntime,nsweep,size_mpi
 153  integer :: rank_to_send
 154  real(dp) :: cutoff,factor,time,tolMSDF,tolMSDS,tolMSDE,tolMSDFS
 155  real(dp),parameter :: HaBohr_meVAng = 27.21138386 / 0.529177249
 156  logical :: iam_master,need_verbose,need_positive,converge
 157  logical :: need_anharmstr,need_spcoupling,ditributed_coefficients
 158  logical :: need_only_odd_power,need_only_even_power,need_initialize_data
 159 !arrays
 160  real(dp) :: mingf(4)
 161  integer :: sc_size(3)
 162  integer,allocatable  :: buffsize(:),buffdisp(:),buffin(:)
 163  integer,allocatable  :: list_coeffs(:),list_coeffs_tmp(:),list_coeffs_tmp2(:)
 164  integer,allocatable  :: my_coeffindexes(:),singular_coeffs(:)
 165  integer,allocatable  :: my_coefflist(:) ,stat_coeff(:)
 166  real(dp),allocatable :: buffGF(:,:),coeff_values(:),energy_coeffs(:,:)
 167  real(dp),allocatable :: energy_coeffs_tmp(:,:)
 168  real(dp),allocatable :: fcart_coeffs(:,:,:,:),gf_values(:,:),gf_mpi(:,:)
 169  real(dp),allocatable :: fcart_coeffs_tmp(:,:,:,:),strten_coeffs_tmp(:,:,:)
 170  real(dp),allocatable :: strten_coeffs(:,:,:)
 171  type(polynomial_coeff_type),allocatable :: my_coeffs(:)
 172  type(polynomial_coeff_type),target,allocatable :: coeffs_tmp(:)
 173  type(polynomial_coeff_type),pointer :: coeffs_in(:)
 174  type(fit_data_type) :: fit_data
 175  character(len=1000) :: message
 176  character(len=fnlen) :: filename
 177  character(len=3)  :: i_char
 178  character(len=7)  :: j_char
 179 ! *************************************************************************
 180 
 181 !MPI variables
 182  master = 0
 183  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
 184  iam_master = (my_rank == master)
 185 
 186 !Initialisation of optional arguments
 187  need_verbose = .TRUE.
 188  if(present(verbose)) need_verbose = verbose
 189  need_initialize_data = .TRUE.
 190  if(present(initialize_data)) need_initialize_data = initialize_data
 191  need_positive = .FALSE.
 192  if(present(positive)) need_positive = positive
 193  need_anharmstr = .FALSE.
 194  if(present(anharmstr)) need_anharmstr = anharmstr
 195  need_spcoupling = .TRUE.
 196  if(present(spcoupling)) need_spcoupling = spcoupling
 197  need_only_odd_power = .FALSE.
 198  if(present(only_odd_power)) need_only_odd_power = only_odd_power
 199  need_only_even_power = .FALSE.
 200  if(present(only_even_power)) need_only_even_power = only_even_power
 201  if(need_only_odd_power.and.need_only_even_power)then
 202       write(message, '(3a)' )&
 203 &       'need_only_odd_power and need_only_even_power are both true',ch10,&
 204 &       'Action: contact abinit group'
 205    MSG_ERROR(message)
 206  end if
 207  max_power_strain_in = 1
 208  if(present(max_power_strain))then
 209    max_power_strain_in = max_power_strain
 210  end if
 211  if(max_power_strain_in <= 0)then
 212       write(message, '(3a)' )&
 213 &       'max_power_strain can not be inferior or equal to zero',ch10,&
 214 &       'Action: contact abinit group'
 215    MSG_ERROR(message)
 216  end if
 217 
 218 !Set the tolerance for the fit
 219  tolMSDF=zero;tolMSDS=zero;tolMSDE=zero;tolMSDFS=zero
 220  if(present(fit_tolMSDF)) tolMSDF  = fit_tolMSDF
 221  if(present(fit_tolMSDS)) tolMSDS  = fit_tolMSDS
 222  if(present(fit_tolMSDE)) tolMSDE  = fit_tolMSDE
 223  if(present(fit_tolMSDFS))tolMSDFS = fit_tolMSDFS
 224 
 225  if(need_verbose) then
 226    write(message,'(a,(80a))') ch10,('=',ii=1,80)
 227    call wrtout(ab_out,message,'COLL')
 228    call wrtout(std_out,message,'COLL')
 229    write(message,'(2a)') ch10,' Starting Fit process'
 230    call wrtout(ab_out,message,'COLL')
 231    call wrtout(std_out,message,'COLL')
 232    write(message,'(a,(80a))') ch10,('-',ii=1,80)
 233    call wrtout(ab_out,message,'COLL')
 234    call wrtout(std_out,message,'COLL')
 235  end if
 236 
 237  ditributed_coefficients = .true.
 238  if(option==2) ditributed_coefficients = .false.
 239 
 240 !if the number of atoms in reference supercell into effpot is not correct,
 241 !wrt to the number of atom in the hist, we set map the hist and set the good supercell
 242  if (size(hist%xred,2) /= eff_pot%supercell%natom) then
 243    call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose)
 244  end if
 245 
 246 !Set the cut off
 247  cutoff = zero
 248  if(present(cutoff_in))then
 249    cutoff = cutoff_in
 250  end if
 251 !If the cutoff is set to zero, we define a default value
 252  if(abs(cutoff)<tol16)then
 253    do ii=1,3
 254      cutoff = cutoff + sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
 255 &                           eff_pot%supercell%rprimd(ii,2)**2+&
 256 &                           eff_pot%supercell%rprimd(ii,3)**2)
 257    end do
 258    cutoff = cutoff / 3.0_dp
 259  end if
 260 !we get the size of the supercell in the hist file
 261  do ii=1,3
 262    sc_size(ii) = int(anint(sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
 263 &                               eff_pot%supercell%rprimd(ii,2)**2+&
 264 &                               eff_pot%supercell%rprimd(ii,3)**2) / &
 265 &                          sqrt(eff_pot%crystal%rprimd(ii,1)**2+&
 266 &                               eff_pot%crystal%rprimd(ii,2)**2+&
 267 &                               eff_pot%crystal%rprimd(ii,3)**2)))
 268  end do
 269 
 270 
 271 !Get the list of coefficients to fit:
 272 !get from the eff_pot type (from the input)
 273 !or
 274 !regenerate the list
 275  my_ncoeff = 0
 276  ncoeff_tot = 0
 277  ncoeff_model = eff_pot%anharmonics_terms%ncoeff
 278 
 279 !Reset ncoeff_tot
 280  if(ncoeff_model > 0)then
 281    if(need_verbose)then
 282      write(message, '(4a)' )ch10,' The coefficients present in the effective',&
 283 &    ' potential will be used for the fit'
 284      call wrtout(std_out,message,'COLL')
 285      call wrtout(ab_out,message,'COLL')
 286    end if
 287  end if
 288 
 289  if(generateterm == 1)then
 290 ! we need to regerate them
 291    if(need_verbose)then
 292      write(message, '(4a)' )ch10,' The coefficients for the fit will be generated'
 293      call wrtout(std_out,message,'COLL')
 294      call wrtout(ab_out,message,'COLL')
 295 
 296      write(message,'(a,F6.3,a)') " Cut-off of ",cutoff," Angstrom is imposed"
 297      call wrtout(std_out,message,'COLL')
 298    end if
 299 
 300    call polynomial_coeff_getNorder(coeffs_tmp,eff_pot%crystal,cutoff,my_ncoeff,ncoeff_tot,power_disps,&
 301 &                                  max_power_strain_in,0,sc_size,comm,anharmstr=need_anharmstr,&
 302 &                                  spcoupling=need_spcoupling,distributed=.true.,&
 303 &                                  only_odd_power=need_only_odd_power,&
 304 &                                  only_even_power=need_only_even_power)
 305  end if
 306 
 307 !Copy the initial coefficients from the model on the CPU 0
 308  ncoeff_tot = ncoeff_tot + ncoeff_model
 309  if(iam_master .and. ncoeff_model > 0) my_ncoeff = my_ncoeff + ncoeff_model
 310 
 311 !Get the list with the number of coeff on each CPU
 312 !In order to be abble to compute the my_coeffindexes array which is for example:
 313 ! if CPU0 has 200  Coeff and CPU1 has 203 Coeff then
 314 ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403
 315 !Also fill the my_coeffs array with the generated coefficients and/or the coefficient from the input xml
 316  ABI_ALLOCATE(buffin,(nproc))
 317  buffin = 0
 318  buffin(my_rank+1) = my_ncoeff
 319  call xmpi_sum(buffin,comm,ierr)
 320  ABI_ALLOCATE(my_coeffindexes,(my_ncoeff))
 321  ABI_ALLOCATE(my_coefflist,(my_ncoeff))
 322  ABI_DATATYPE_ALLOCATE(my_coeffs,(my_ncoeff))
 323  do icoeff=1,my_ncoeff
 324 
 325    jcoeff = icoeff
 326    my_coefflist(icoeff) = icoeff
 327 
 328    if(my_rank==0) then
 329      my_coeffindexes(icoeff) = icoeff
 330    else
 331      my_coeffindexes(icoeff) = sum(buffin(1:my_rank)) + icoeff
 332    end if
 333 
 334 !  Only copy the input coefficients on the CPU0
 335    if(my_rank==0) then
 336      if(icoeff <= ncoeff_model)then
 337        coeffs_in => eff_pot%anharmonics_terms%coefficients
 338      else
 339        coeffs_in => coeffs_tmp
 340        jcoeff = jcoeff-ncoeff_model
 341      end if
 342    else
 343      coeffs_in => coeffs_tmp
 344    end if
 345    call polynomial_coeff_init(one,coeffs_in(jcoeff)%nterm,&
 346 &                             my_coeffs(icoeff),coeffs_in(jcoeff)%terms,&
 347 &                             coeffs_in(jcoeff)%name,&
 348 &                             check=.true.)
 349    call polynomial_coeff_free(coeffs_in(jcoeff))
 350  end do
 351 
 352 !Deallocation
 353  if(allocated(coeffs_tmp)) then
 354    ABI_DATATYPE_DEALLOCATE(coeffs_tmp)
 355  end if
 356  NULLIFY(coeffs_in)
 357  ABI_DEALLOCATE(buffin)
 358 
 359  !wait everybody
 360  call xmpi_barrier(comm)
 361 
 362 !Write the XML with the coefficient before the fit process
 363  if(iam_master)then
 364    filename = "terms_set.xml"
 365 !   call polynomial_coeff_writeXML(my_coeffs,my_ncoeff,filename=filename,newfile=.true.)
 366  end if
 367 
 368 !Reset the output (we free the memory)
 369  call effective_potential_freeCoeffs(eff_pot)
 370 
 371 
 372 !Check if ncycle_in is not zero or superior to ncoeff_tot
 373  if(need_verbose.and.(ncycle_in > ncoeff_tot).or.(ncycle_in<0.and.nfixcoeff /= -1)) then
 374    write(message, '(6a,I0,3a)' )ch10,&
 375 &        ' --- !WARNING',ch10,&
 376 &        '     The number of cycle requested in the input is not correct.',ch10,&
 377 &        '     This number will be set to the maximum of coefficients: ',ncoeff_tot,ch10,&
 378 &        ' ---',ch10
 379      call wrtout(std_out,message,"COLL")
 380    end if
 381 
 382 !Use fixcoeff
 383 !ncycle_tot store the curent number of coefficient in the model
 384 !Do not reset this variable...
 385  ncycle_tot = 0
 386  if (nfixcoeff == -1)then
 387    write(message, '(3a)')' nfixcoeff is set to -1, the coefficients present in the model',&
 388 &                        ' are imposed.',ch10
 389    ncycle_tot = ncycle_tot + ncoeff_model
 390  else
 391    if (nfixcoeff > 0)then
 392      if(maxval(fixcoeff(:)) > ncoeff_tot) then
 393        write(message, '(4a,I0,6a)' )ch10,&
 394 &        ' --- !WARNING',ch10,&
 395 &        '     The value ',maxval(fixcoeff(:)),' is not in the list.',ch10,&
 396 &        '     Start from scratch...',ch10,&
 397 &        ' ---',ch10
 398      else
 399        ncycle_tot = ncycle_tot + nfixcoeff
 400        write(message, '(2a)')' Some coefficients are imposed from the input.',ch10
 401      end if
 402    else
 403      write(message, '(4a)')' There is no coefficient imposed from the input.',ch10,&
 404 &                        ' Start from scratch',ch10
 405    end if
 406  end if
 407 
 408  if(need_verbose) call wrtout(std_out,message,'COLL')
 409 
 410 !Compute the number of cycle:
 411  ncycle     = ncycle_in
 412 !Compute the maximum number of cycle
 413  ncycle_max = ncycle_in + ncycle_tot
 414 
 415 !Check if the number of request cycle + the initial number of coeff is superior to
 416  !the maximum number of coefficient allowed
 417  if(ncycle_max > ncoeff_tot) then
 418    ncycle = ncoeff_tot - ncycle_tot
 419    ncycle_max = ncoeff_tot
 420    write(message, '(4a,I0,2a,I0,2a,I0,3a)' )ch10,&
 421 &      ' --- !WARNING',ch10,&
 422 &      '     The number of cycle + the number of imposed coefficients: ',ncycle_max,ch10,&
 423 &      '     is superior to the maximum number of coefficients in the initial list: ',ncoeff_tot,ch10,&
 424 &      '     The number of cycle is set to ',ncycle,ch10,&
 425 &      ' ---',ch10
 426    if(need_verbose) call wrtout(std_out,message,'COLL')
 427  else if (option==2)then
 428 !  Always set to the maximum
 429    ncycle_max = ncoeff_tot
 430  end if
 431 
 432 !Initialisation of constants
 433  ntime    = hist%mxhist
 434  natom_sc = eff_pot%supercell%natom
 435  ncell    = eff_pot%supercell%ncells
 436  factor   = 1._dp/natom_sc
 437 
 438 !Initialisation of arrays:
 439  ABI_ALLOCATE(energy_coeffs_tmp,(ncycle_max,ntime))
 440  ABI_ALLOCATE(list_coeffs,(ncycle_max))
 441  ABI_ALLOCATE(fcart_coeffs_tmp,(3,natom_sc,ncycle_max,ntime))
 442  ABI_ALLOCATE(strten_coeffs_tmp,(6,ntime,ncycle_max))
 443  list_coeffs  = 0
 444 
 445 !if ncycle_tot > 0 fill list_coeffs with the fixed coefficients
 446  if(ncycle_tot > 0)then
 447    do ii = 1,ncycle_tot
 448      if(nfixcoeff == -1)then
 449        if(ii <= ncoeff_model)then
 450          list_coeffs(ii) = ii
 451        end if
 452      else
 453        list_coeffs(ii) = fixcoeff(ii)
 454      end if
 455    end do
 456  end if
 457 
 458 !Get the decomposition for each coefficients of the forces and stresses for
 459 !each atoms and each step  equations 11 & 12 of  PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
 460  if(need_verbose)then
 461    write(message, '(a)' ) ' Initialisation of the fit process...'
 462    call wrtout(std_out,message,'COLL')
 463  end if
 464 !Before the fit, compute constants with fit_data_compute.
 465 !Conpute the strain of each configuration.
 466 !Compute the displacmeent of each configuration.
 467 !Compute the variation of the displacement due to strain of each configuration.
 468 !Compute fixed forces and stresse and get the standard deviation.
 469 !Compute Sheppard and al Factors  \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
 470  call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose)
 471 
 472 !Get the decomposition for each coefficients of the forces,stresses and energy for
 473 !each atoms and each step  (see equations 11 & 12 of  
 474 !PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]]+ allocation
 475 !If the user does not turn off this initialization, we store all the informations for the fit,
 476 !it will reduce the computation time but increase a lot the memory...
 477  if(need_initialize_data)then
 478    ABI_ALLOCATE(energy_coeffs,(my_ncoeff,ntime))
 479    ABI_ALLOCATE(fcart_coeffs,(3,natom_sc,my_ncoeff,ntime))
 480    ABI_ALLOCATE(strten_coeffs,(6,ntime,my_ncoeff))
 481    call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
 482 &                                 fit_data%training_set%displacement,&
 483 &                                 energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
 484 &                                 my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
 485 &                                 strten_coeffs,fit_data%training_set%ucvol,my_coefflist,my_ncoeff)
 486  else
 487 !  Allocate just 1 dimension ! Save MEMORY !
 488    ABI_ALLOCATE(energy_coeffs,(1,ntime))
 489    ABI_ALLOCATE(fcart_coeffs,(3,natom_sc,1,ntime))
 490    ABI_ALLOCATE(strten_coeffs,(6,ntime,1))   
 491  end if
 492 
 493 !Allocation of arrays
 494  ABI_DATATYPE_ALLOCATE(coeffs_tmp,(ncycle_max))
 495  ABI_ALLOCATE(singular_coeffs,(max(1,my_ncoeff)))
 496  ABI_ALLOCATE(coeff_values,(ncycle_max))
 497  ABI_ALLOCATE(gf_values,(4,max(1,my_ncoeff)))
 498  ABI_ALLOCATE(list_coeffs_tmp,(ncycle_max))
 499  ABI_ALLOCATE(list_coeffs_tmp2,(ncycle_max))
 500  ABI_ALLOCATE(stat_coeff,(ncoeff_tot))
 501  coeff_values = zero
 502  singular_coeffs = 0
 503  stat_coeff = 0
 504 !Set mpi buffer
 505 !Set the bufsize for mpi allgather
 506  ABI_ALLOCATE(buffsize,(nproc))
 507  ABI_ALLOCATE(buffdisp,(nproc))
 508  ABI_ALLOCATE(buffGF,(5,1))
 509  ABI_ALLOCATE(gf_mpi,(5,nproc))
 510  buffsize(:) = 0
 511  buffdisp(1) = 0
 512  do ii= 1,nproc
 513    buffsize(ii) =  5
 514  end do
 515  do ii = 2,nproc
 516    buffdisp(ii) = buffdisp(ii-1) + buffsize(ii-1)
 517  end do
 518  size_mpi = 5*nproc
 519 !If some coeff are imposed by the input, we need to fill the arrays
 520 !with this coeffs and broadcast to the others CPUs :
 521  if(ncycle_tot>=1)then
 522    do icycle = 1,ncycle_tot
 523      list_coeffs_tmp(icycle) = icycle
 524      rank_to_send = 0
 525      do icoeff=1,my_ncoeff
 526        if((my_coeffindexes(icoeff)==list_coeffs(icycle)))then
 527 
 528          if(need_initialize_data)then
 529            my_icoeff = icoeff
 530          else
 531            my_icoeff = 1
 532 !          Need to initialized the data for the fit for this coefficient 
 533            call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
 534 &                                          fit_data%training_set%displacement,&
 535 &                                          energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
 536 &                                          my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
 537 &                                          strten_coeffs,fit_data%training_set%ucvol,&
 538 &                                          my_coefflist(icoeff),1)
 539          end if
 540          
 541          energy_coeffs_tmp(icycle,:)    = energy_coeffs(my_icoeff,:)
 542          fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:)
 543          strten_coeffs_tmp(:,:,icycle)  = strten_coeffs(:,:,my_icoeff)
 544          rank_to_send = my_rank
 545          call polynomial_coeff_free(coeffs_tmp(icycle))
 546          call polynomial_coeff_init(coeff_values(icycle),my_coeffs(icoeff)%nterm,&
 547 &                                   coeffs_tmp(icycle),my_coeffs(icoeff)%terms,&
 548 &                                   my_coeffs(icoeff)%name,&
 549 &                                   check=.false.)
 550          exit
 551        end if
 552      end do
 553 !    Need to send the rank with the chosen coefficient
 554      call xmpi_sum(rank_to_send, comm, ierr)
 555 !    Boadcast the coefficient
 556      call xmpi_bcast(energy_coeffs_tmp(icycle,:), rank_to_send, comm, ierr)
 557      call xmpi_bcast(fcart_coeffs_tmp(:,:,icycle,:) , rank_to_send, comm, ierr)
 558      call xmpi_bcast(strten_coeffs_tmp(:,:,icycle), rank_to_send, comm, ierr)
 559      call polynomial_coeff_broadcast(coeffs_tmp(icycle), rank_to_send, comm)
 560    end do
 561  end if
 562 
 563 !Waiting for all
 564  if(nproc > 1)  then
 565    if(need_verbose)then
 566      write(message, '(a)') ' Initialisation done... waiting for all the CPU'
 567      call wrtout(std_out,message,'COLL')
 568    end if
 569    call xmpi_barrier(comm)
 570  end if
 571 
 572 !Compute GF, coeff_values,strten_coeffs and fcart_coeffs are set to zero
 573 !it means that only the harmonic part wiil be computed
 574  coeff_values = zero
 575 
 576  call fit_polynomial_coeff_computeGF(coeff_values,energy_coeffs,fit_data%energy_diff,fcart_coeffs,&
 577 &                                    fit_data%fcart_diff,gf_values(:,1),int((/1/)),natom_sc,&
 578 &                                    0,my_ncoeff,ntime,strten_coeffs,fit_data%strten_diff,&
 579 &                                    fit_data%training_set%sqomega)
 580 
 581 !Print the standard deviation before the fit
 582  write(message,'(3a,ES24.16,4a,ES24.16,2a,ES24.16,2a,ES24.16,a)' ) &
 583 &                   ' Mean Standard Deviation values at the begining of the fit process (meV/atm):',&
 584 &               ch10,'   Energy          : ',&
 585 &               gf_values(4,1)*Ha_EV*1000*factor  ,ch10,&
 586 &                    ' Goal function values at the begining of the fit process (eV^2/A^2):',ch10,&
 587 &                    '   Forces+Stresses : ',&
 588 &               gf_values(1,1)*(HaBohr_meVAng)**2,ch10,&
 589 &                    '   Forces          : ',&
 590 &               gf_values(2,1)*(HaBohr_meVAng)**2,ch10,&
 591 &                    '   Stresses        : ',&
 592 &               gf_values(3,1)*(HaBohr_meVAng)**2,ch10
 593  if(need_verbose)then
 594    call wrtout(ab_out,message,'COLL')
 595    call wrtout(std_out,message,'COLL')
 596  end if
 597 
 598  select case(option)
 599 
 600  case(1)
 601    !Option 1, we select the coefficients one by one
 602    if(need_verbose.and.ncycle > 0)then
 603      write(message,'(a,3x,a,10x,a,14x,a,14x,a,14x,a)') " N","Selecting","MSDE","MSDFS","MSDF","MSDS"
 604      call wrtout(ab_out,message,'COLL')
 605      write(message,'(4x,a,6x,a,8x,a,8x,a,8x,a)') "Coefficient","(meV/atm)","(eV^2/A^2)","(eV^2/A^2)",&
 606 &                                            "(eV^2/A^2)"
 607      call wrtout(ab_out,message,'COLL')
 608    end if
 609 
 610 !  Start fit process
 611    do icycle_tmp = 1,ncycle
 612      icycle = ncycle_tot + 1
 613      list_coeffs_tmp(icycle) = icycle
 614      if(need_verbose)then
 615        write(message, '(4a,I0,a)')ch10,'--',ch10,' Try to find the best model with ',&
 616 &                                 icycle,' coefficient'
 617        if(icycle > 1)  write(message, '(2a)') trim(message),'s'
 618        if(nproc > 1)  then
 619          if(my_ncoeff>=1) then
 620            write(message, '(2a,I0,a)')trim(message), ' (only the ',my_ncoeff,&
 621 &                                                ' first are printed for this CPU)'
 622          else
 623            write(message, '(2a)')trim(message), ' (no coefficient treated by this CPU)'
 624          end if
 625        end if
 626        call wrtout(std_out,message,'COLL')
 627        if(icycle>1 .or. any(list_coeffs(:) > zero))then
 628          write(message, '(3a)') ' The coefficient numbers from the previous cycle are:',ch10,' ['
 629          do ii=1,icycle-1
 630            if(ii<icycle-1)then
 631              write(message, '(a,I0,a)') trim(message),list_coeffs(ii),','
 632            else
 633              write(message, '(a,I0)') trim(message),list_coeffs(ii)
 634            end if
 635          end do
 636          write(message, '(3a)') trim(message),']',ch10
 637          call wrtout(std_out,message,'COLL')
 638        end if
 639 
 640        write(message,'(2x,a,12x,a,14x,a,13x,a,14x,a)') " Testing","MSDE","MSDFS","MSDF","MSDS"
 641        call wrtout(std_out,message,'COLL')
 642        write(message,'(a,7x,a,8x,a,8x,a,8x,a)') " Coefficient","(meV/atm)","(eV^2/A^2)","(eV^2/A^2)",&
 643 &                                            "(eV^2/A^2)"
 644        call wrtout(std_out,message,'COLL')
 645      end if!End if verbose
 646 
 647 !    Reset gf_values
 648      gf_values(:,:) = zero
 649      do icoeff=1,my_ncoeff
 650 !    cycle if this coefficient is not allowed
 651        if(any(list_coeffs==my_coeffindexes(icoeff)).or.singular_coeffs(icoeff) == 1) cycle
 652        if(nbancoeff >= 1)then
 653          if(any(bancoeff==my_coeffindexes(icoeff))) cycle
 654        end if
 655        list_coeffs(icycle) = my_coeffindexes(icoeff)
 656 
 657        if(need_initialize_data)then
 658          my_icoeff = icoeff
 659        else
 660 !        Need to initialized the data for the fit for this coefficient
 661          my_icoeff = 1        
 662          call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
 663 &                                        fit_data%training_set%displacement,&
 664 &                                        energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
 665 &                                        my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
 666 &                                        strten_coeffs,fit_data%training_set%ucvol,&
 667 &                                        my_coefflist(icoeff),1)
 668        end if
 669        
 670 !      Fill the temporary arrays
 671        energy_coeffs_tmp(icycle,:)    = energy_coeffs(my_icoeff,:)
 672        fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:)
 673        strten_coeffs_tmp(:,:,icycle)  = strten_coeffs(:,:,my_icoeff)
 674 
 675 !      call the fit process routine
 676 !      This routine solves the linear system proposed 
 677 !      by C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
 678        call fit_polynomial_coeff_solve(coeff_values(1:icycle),fcart_coeffs_tmp,fit_data%fcart_diff,&
 679 &                                      energy_coeffs_tmp,fit_data%energy_diff,info,&
 680 &                                      list_coeffs_tmp(1:icycle),natom_sc,icycle,ncycle_max,ntime,&
 681 &                                      strten_coeffs_tmp,fit_data%strten_diff,&
 682 &                                      fit_data%training_set%sqomega)
 683        if(info==0)then
 684          if (need_positive.and.any(coeff_values(nfixcoeff+1:icycle) < zero)) then
 685            write(message, '(a)') ' Negative value detected...'
 686            gf_values(:,icoeff) = zero
 687            coeff_values = zero
 688          else
 689            call fit_polynomial_coeff_computeGF(coeff_values(1:icycle),energy_coeffs_tmp,&
 690 &                                            fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
 691 &                                            gf_values(:,icoeff),list_coeffs_tmp(1:icycle),natom_sc,&
 692 &                                            icycle,ncycle_max,ntime,strten_coeffs_tmp,&
 693 &                                            fit_data%strten_diff,fit_data%training_set%sqomega)
 694 
 695 
 696            write (j_char, '(i7)') my_coeffindexes(icoeff)
 697            write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char),&
 698 &                                   gf_values(4,icoeff)* 1000*Ha_ev *factor,&
 699 &                                   gf_values(1,icoeff)*HaBohr_meVAng**2,&
 700 &                                   gf_values(2,icoeff)*HaBohr_meVAng**2,&
 701 &                                   gf_values(3,icoeff)*HaBohr_meVAng**2
 702          end if
 703        else!In this case the matrix is singular
 704          gf_values(:,icoeff) = zero
 705          singular_coeffs(icoeff) = 1
 706          write(message, '(a)') ' The matrix is singular...'
 707        end if
 708        if(need_verbose) call wrtout(std_out,message,'COLL')
 709      end do
 710 
 711 !    find the best coeff on each CPU
 712      mingf(:)  = 9D99
 713      index_min = 0
 714      do icoeff=1,my_ncoeff
 715        if(gf_values(1,icoeff) < zero) cycle
 716        if(abs(gf_values(1,icoeff)) <tol16) cycle
 717        if(gf_values(1,icoeff) < mingf(1) ) then
 718          mingf(:) = gf_values(:,icoeff)
 719          index_min = my_coeffindexes(icoeff)
 720        end if
 721      end do
 722 
 723 !    MPI GATHER THE BEST COEFF ON EACH CPU
 724      if(nproc > 1)then
 725        buffGF(1,1) = index_min
 726        buffGF(2:5,1) =  mingf(:)
 727 
 728        call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr)
 729 !      find the best coeff
 730        mingf(:)    = 9D99
 731        index_min= 0
 732        do icoeff=1,nproc
 733          if(gf_mpi(2,icoeff) < zero) cycle
 734          if(abs(gf_mpi(2,icoeff)) < tol16) cycle
 735          if(gf_mpi(2,icoeff) < mingf(1) ) then
 736            mingf(:) = gf_mpi(2:5,icoeff)
 737            index_min = int(gf_mpi(1,icoeff))
 738          end if
 739        end do
 740      end if
 741 
 742 !    Check if there is still coefficient
 743      if(index_min==0) then
 744        exit
 745      else
 746        list_coeffs(icycle) = index_min
 747      end if
 748 
 749 !    Check if this coeff is treat by this cpu and fill the
 750 !    temporary array before broadcast
 751      rank_to_send = 0
 752      do icoeff=1,my_ncoeff
 753 
 754 
 755        if((my_coeffindexes(icoeff)==list_coeffs(icycle)))then
 756 
 757          if(need_initialize_data)then
 758            my_icoeff = icoeff
 759          else
 760 !          Need to initialized the data for the fit for this coefficient
 761            my_icoeff = 1           
 762            call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
 763 &                                          fit_data%training_set%displacement,&
 764 &                                          energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
 765 &                                          my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
 766 &                                          strten_coeffs,fit_data%training_set%ucvol,&
 767 &                                          my_coefflist(icoeff),1)
 768          end if
 769 
 770          energy_coeffs_tmp(icycle,:)    = energy_coeffs(my_icoeff,:)
 771          fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:)
 772          strten_coeffs_tmp(:,:,icycle)  = strten_coeffs(:,:,my_icoeff)
 773          call polynomial_coeff_free(coeffs_tmp(icycle))
 774          call polynomial_coeff_init(coeff_values(icycle),my_coeffs(icoeff)%nterm,&
 775 &                                   coeffs_tmp(icycle),my_coeffs(icoeff)%terms,&
 776 &                                   my_coeffs(icoeff)%name,&
 777 &                                   check=.false.)
 778 
 779          rank_to_send = my_rank
 780          exit
 781        end if
 782      end do
 783 !    Need to send the rank with the chosen coefficient
 784      call xmpi_sum(rank_to_send, comm, ierr)
 785 
 786 !    Boadcast the coefficient
 787      call xmpi_bcast(energy_coeffs_tmp(icycle,:), rank_to_send, comm, ierr)
 788      call xmpi_bcast(fcart_coeffs_tmp(:,:,icycle,:) , rank_to_send, comm, ierr)
 789      call xmpi_bcast(strten_coeffs_tmp(:,:,icycle), rank_to_send, comm, ierr)
 790      call polynomial_coeff_broadcast(coeffs_tmp(icycle), rank_to_send, comm)
 791 
 792      if(need_verbose) then
 793        write(message, '(a,I0,2a)' )' Selecting the coefficient number ',list_coeffs(icycle),&
 794 &                                   ' ===> ',trim(coeffs_tmp(icycle)%name)
 795        call wrtout(std_out,message,'COLL')
 796 
 797        write(message, '(2a,I0,a,ES24.16)' )' Standard deviation of the energy for',&
 798 &                                        ' the iteration ',icycle_tmp,' (meV/atm): ',&
 799 &                         mingf(4)* Ha_eV *1000 *factor
 800        call wrtout(std_out,message,'COLL')
 801 
 802        write (i_char, '(i3)') icycle
 803        write (j_char, '(i7)') list_coeffs(icycle)
 804        write(message, '(a,a,3x,a,3x,4ES18.10)') " ",adjustl(i_char),adjustl(j_char),&
 805 &                                    mingf(4)* 1000*Ha_eV *factor,&
 806 &                                    mingf(1)*HaBohr_meVAng**2,&
 807 &                                    mingf(2)*HaBohr_meVAng**2,&
 808 &                                    mingf(3)*HaBohr_meVAng**2
 809        call wrtout(ab_out,message,'COLL')
 810      end if
 811 
 812      ncycle_tot = ncycle_tot + 1
 813 
 814 !    Check the stopping criterion
 815      converge = .false.
 816      if(tolMSDE  > zero)then
 817        if(abs(tolMSDE) > abs(mingf(4)*1000*Ha_eV *factor))then
 818          write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
 819 &                                                mingf(4)*1000*Ha_eV * factor ," < ",tolMSDE,&
 820 &                                              ' for MSDE'
 821          converge = .true.
 822        end if
 823      end if
 824      if(tolMSDF  > zero) then
 825        if(abs(tolMSDF) > abs(mingf(2)*HaBohr_meVAng**2))then
 826          write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
 827 &                                                  mingf(2)*HaBohr_meVAng**2 ," < ",tolMSDF,&
 828 &                                              ' for MSDF'
 829          converge = .true.
 830        end if
 831      end if
 832      if(tolMSDS  > zero) then
 833        if(abs(tolMSDS) > abs(mingf(3)*HaBohr_meVAng**2))then
 834          write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
 835 &                                                  mingf(3)*HaBohr_meVAng**2 ," < ",tolMSDS,&
 836 &                                              ' for MSDS'
 837          converge = .true.
 838        end if
 839      end if
 840      if(tolMSDFS > zero)then
 841        if(abs(tolMSDFS) > abs(mingf(1)*HaBohr_meVAng**2))then
 842          write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
 843 &                                                  mingf(1)*HaBohr_meVAng**2 ," < ",tolMSDFS,&
 844 &                                              ' for MSDFS'
 845          converge = .true.
 846        end if
 847      end if
 848      if(converge)then
 849        call wrtout(ab_out,message,'COLL')
 850        call wrtout(std_out,message,'COLL')
 851        exit
 852      else
 853        if(any((/abs(tolMSDE),abs(tolMSDF),abs(tolMSDS),abs(tolMSDFS)/) > tol20) .and.&
 854 &         icycle_tmp == ncycle)then
 855          write(message,'(2a,I0,a)') ch10," WARNING: ",ncycle,&
 856 &                                   " cycles was not enougth to converge the fit process"
 857          call wrtout(ab_out,message,'COLL')
 858          call wrtout(std_out,message,'COLL')
 859        end if
 860      end if
 861 
 862    end do
 863 
 864  case(2)
 865 
 866 !  Monte Carlo selection
 867    nsweep = 10000
 868 !  If no coefficient imposed in the inputs we reset the goal function
 869    if (ncycle_tot == 0) then
 870      gf_values(:,:) = zero
 871      mingf(:) = 9D99
 872    else
 873      mingf = gf_values(:,1)
 874    end if
 875    call cpu_time(time)
 876    call ZBQLINI(int(time*1000000/(my_rank+1)))
 877    if(need_verbose)then
 878      write(message,'(a,I0,a)') " Start Monte Carlo simulations on ", nproc," CPU"
 879      if(nproc>1) write(message,'(2a)') trim(message)," (only print result of the master)"
 880      call wrtout(std_out,message,'COLL')
 881      call wrtout(ab_out,message,'COLL')
 882      write(message,'(a,2x,a,9x,a,14x,a,13x,a,14x,a)') ch10," Iteration ","MSDE","MSDFS","MSDF","MSdS"
 883      call wrtout(std_out,message,'COLL')
 884      write(message,'(a,5x,a,8x,a,8x,a,8x,a)') "              ","(meV/atm)","(eV^2/A^2)","(eV^2/A^2)",&
 885 &                                            "(eV^2/A^2)"
 886      call wrtout(std_out,message,'COLL')
 887 
 888    end if
 889 
 890    do ii = 1,1!nyccle
 891      do isweep =1,nsweep
 892        write (j_char, '(i7)') isweep
 893 !TEST_AM
 894        icycle_tmp = int(ZBQLU01(zero)*(ncycle+1-1))+1
 895        icycle_tmp = ncycle
 896        do icycle=1,icycle_tmp
 897          icoeff = int(ZBQLU01(zero)*(my_ncoeff))+1
 898 !         icycle = int(ZBQLU01(zero)*(ncycle))+1
 899          list_coeffs_tmp2(icycle) = icoeff
 900          list_coeffs_tmp(icycle)= icycle
 901 !        Fill the temporary arrays
 902          energy_coeffs_tmp(icycle,:)    = energy_coeffs(icoeff,:)
 903          fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,icoeff,:)
 904          strten_coeffs_tmp(:,:,icycle)  = strten_coeffs(:,:,icoeff)
 905        end do
 906 !TEST_AM
 907 
 908 !      call the fit process routine
 909 !      This routine solves the linear system proposed by 
 910 !      C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
 911        call fit_polynomial_coeff_solve(coeff_values(1:icycle_tmp),fcart_coeffs_tmp,fit_data%fcart_diff,&
 912 &                                      energy_coeffs_tmp,fit_data%energy_diff,info,&
 913 &                                      list_coeffs_tmp(1:icycle_tmp),natom_sc,icycle_tmp,ncycle_max,&
 914 &                                      ntime,strten_coeffs_tmp,fit_data%strten_diff,&
 915 &                                      fit_data%training_set%sqomega)
 916        if(info==0)then
 917          call fit_polynomial_coeff_computeGF(coeff_values(1:icycle_tmp),energy_coeffs_tmp,&
 918 &                                            fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
 919 &                                            gf_values(:,1),list_coeffs_tmp(1:icycle_tmp),natom_sc,&
 920 &                                            icycle_tmp,ncycle_max,ntime,strten_coeffs_tmp,&
 921 &                                            fit_data%strten_diff,fit_data%training_set%sqomega)
 922 
 923        else!In this case the matrix is singular
 924          gf_values(:,icoeff) = zero
 925          singular_coeffs(icoeff) = 1
 926        end if
 927 
 928        if(gf_values(1,1) > zero.and.abs(gf_values(1,1))>tol16.and.&
 929 &         gf_values(1,1) < mingf(1) ) then
 930          mingf = gf_values(:,1)
 931          list_coeffs(1:icycle_tmp) = list_coeffs_tmp2(1:icycle_tmp)
 932          ncycle_tot = icycle_tmp
 933 
 934          write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char),&
 935 &                                   gf_values(4,1)* 1000*Ha_ev *factor,&
 936 &                                   gf_values(1,1)*HaBohr_meVAng**2,&
 937 &                                   gf_values(2,1)*HaBohr_meVAng**2,&
 938 &                                   gf_values(3,1)*HaBohr_meVAng**2
 939          if(need_verbose) call wrtout(std_out,message,'COLL')
 940        else
 941          list_coeffs_tmp2(1:icycle_tmp) = list_coeffs(1:icycle_tmp)
 942        end if
 943      end do
 944 
 945      if(nproc > 1) then
 946 !TEST_AM
 947        do iproc=1,ncycle_tot
 948          stat_coeff(list_coeffs(iproc)) =  stat_coeff(list_coeffs(iproc)) + 1
 949        end do
 950 !TEST_AM
 951 
 952 !    Find the best model on all the CPUs
 953        buffGF(1,1) = zero
 954        buffGF(2:5,1) =  mingf(:)
 955        call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr)
 956 !      find the best coeff
 957        mingf(:) = 9D99
 958        index_min= 0
 959        do iproc=1,nproc
 960          if(gf_mpi(2,iproc) < zero) cycle
 961          if(abs(gf_mpi(2,iproc)) <tol16) cycle
 962          if(gf_mpi(2,iproc) < mingf(1) ) then
 963            mingf(:) = gf_mpi(2:5,iproc)
 964            index_min = int(gf_mpi(1,iproc))
 965            rank_to_send = iproc-1
 966          end if
 967        end do
 968        write(message, '(2a,I0)') ch10,' Best model found on the CPU: ', rank_to_send
 969        call wrtout(std_out,message,'COLL')
 970      end if
 971    end do
 972 
 973 !TEST_AM
 974    call xmpi_sum(stat_coeff, comm, ierr)
 975    do ii=1,ncoeff_tot
 976      write(100,*) ii,stat_coeff(ii)
 977    end do
 978    close(100)
 979 !TEST_AM
 980 
 981 !  Transfert final model
 982    if(nproc>1)then
 983      call xmpi_bcast(ncycle_tot,rank_to_send,comm,ierr)
 984      call xmpi_bcast(list_coeffs(1:ncycle_tot),rank_to_send,comm,ierr)
 985    end if
 986    do ii=1,ncycle_tot
 987      icoeff = list_coeffs(ii)
 988      list_coeffs_tmp(ii) = ii
 989 !    Fill the temporary arrays
 990      energy_coeffs_tmp(ii,:)    = energy_coeffs(icoeff,:)
 991      fcart_coeffs_tmp(:,:,ii,:) = fcart_coeffs(:,:,icoeff,:)
 992      strten_coeffs_tmp(:,:,ii)  = strten_coeffs(:,:,icoeff)
 993      call polynomial_coeff_free(coeffs_tmp(ii))
 994      call polynomial_coeff_init(one,my_coeffs(icoeff)%nterm,&
 995 &                               coeffs_tmp(ii),my_coeffs(icoeff)%terms,&
 996 &                               my_coeffs(icoeff)%name,&
 997 &                               check=.false.)
 998    end do
 999  end select
1000 
1001 !This routine solves the linear system proposed by 
1002 ! C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
1003  if(ncycle_tot > 0)then
1004 
1005    call fit_polynomial_coeff_solve(coeff_values(1:ncycle_tot),fcart_coeffs_tmp,fit_data%fcart_diff,&
1006 &                                  energy_coeffs_tmp,fit_data%energy_diff,info,&
1007 &                                  list_coeffs_tmp(1:ncycle_tot),natom_sc,&
1008 &                                  ncycle_tot,ncycle_max,ntime,strten_coeffs_tmp,&
1009 &                                  fit_data%strten_diff,fit_data%training_set%sqomega)
1010 
1011    if(need_verbose) then
1012      write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit process: '
1013      call wrtout(ab_out,message,'COLL')
1014      call wrtout(std_out,message,'COLL')
1015    end if
1016    do ii = 1,ncycle_tot
1017      if(list_coeffs(ii) ==0) cycle
1018 !    Set the value of the coefficient
1019      coeffs_tmp(ii)%coefficient = coeff_values(ii)
1020      if(need_verbose) then
1021        write(message, '(a,I0,a,ES19.10,2a)') " ",list_coeffs(ii)," =>",coeff_values(ii),&
1022          &                                " ",trim(coeffs_tmp(ii)%name)
1023        call wrtout(ab_out,message,'COLL')
1024        call wrtout(std_out,message,'COLL')
1025      end if
1026    end do
1027 
1028    call fit_polynomial_coeff_computeGF(coeff_values(1:ncycle_tot),energy_coeffs_tmp,&
1029 &                                      fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
1030 &                                      gf_values(:,1),list_coeffs_tmp(1:ncycle_tot),natom_sc,&
1031 &                                      ncycle_tot,ncycle_max,ntime,strten_coeffs_tmp,&
1032 &                                      fit_data%strten_diff,fit_data%training_set%sqomega)
1033 
1034    if(need_verbose) then
1035 !  Print the standard deviation after the fit
1036      write(message,'(4a,ES24.16,4a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
1037 &                    ' Mean Standard Deviation values at the end of the fit process (meV/atm):',&
1038 &               ch10,'   Energy          : ',&
1039 &               gf_values(4,1)*Ha_EV*1000*factor ,ch10,&
1040 &                    ' Goal function values at the end of the fit process (eV^2/A^2):',ch10,&
1041 &                    '   Forces+Stresses : ',&
1042 &               gf_values(1,1)*(HaBohr_meVAng)**2,ch10,&
1043 &                    '   Forces          : ',&
1044 &               gf_values(2,1)*(HaBohr_meVAng)**2,ch10,&
1045 &                    '   Stresses        : ',&
1046 &               gf_values(3,1)*(HaBohr_meVAng)**2,ch10
1047      call wrtout(ab_out,message,'COLL')
1048      call wrtout(std_out,message,'COLL')
1049    end if
1050 
1051 !  Set the final set of coefficients into the eff_pot type
1052    call effective_potential_setCoeffs(coeffs_tmp(1:ncycle_tot),eff_pot,ncycle_tot)
1053    call fit_polynomial_coeff_computeMSD(eff_pot,hist,gf_values(4,1),gf_values(2,1),gf_values(1,1),&
1054 &                                       natom_sc,ntime,fit_data%training_set%sqomega,&
1055 &                                       compute_anharmonic=.TRUE.,print_file=.TRUE.)
1056 
1057 !    if(need_verbose) then
1058 ! !  Print the standard deviation after the fit
1059 !      write(message,'(4a,ES24.16,4a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
1060 ! &                    ' Mean Standard Deviation values at the end of the fit process (meV/f.u.):',&
1061 ! &               ch10,'   Energy          : ',&
1062 ! &               gf_values(4,1)*Ha_EV*1000/ ncell ,ch10,&
1063 ! &                    ' Goal function values at the end of the fit process (eV^2/A^2):',ch10,&
1064 ! &                    '   Forces+Stresses : ',&
1065 ! &               (gf_values(1,1)+gf_values(2,1))*(HaBohr_meVAng)**2,ch10,&
1066 ! &                    '   Forces          : ',&
1067 ! &               gf_values(2,1)*(HaBohr_meVAng)**2,ch10,&
1068 ! &                    '   Stresses        : ',&
1069 ! &               gf_values(3,1)*(HaBohr_meVAng)**2,ch10
1070 !      call wrtout(ab_out,message,'COLL')
1071 !      call wrtout(std_out,message,'COLL')
1072 !    end if
1073 
1074  else
1075    if(need_verbose) then
1076      write(message, '(9a)' )ch10,&
1077 &          ' --- !WARNING',ch10,&
1078 &          '     The fit process does not provide possible terms.',ch10,&
1079 &          '     Please make sure than the terms set is correct',ch10,&
1080 &          ' ---',ch10
1081      call wrtout(ab_out,message,'COLL')
1082      call wrtout(std_out,message,'COLL')
1083    end if
1084  end if
1085 !Deallocation of arrays
1086  call fit_data_free(fit_data)
1087 
1088 !Deallocate the temporary coefficient
1089  do ii=1,ncycle_max
1090    call polynomial_coeff_free(coeffs_tmp(ii))
1091  end do
1092  ABI_DATATYPE_DEALLOCATE(coeffs_tmp)
1093  do ii=1,my_ncoeff
1094    call polynomial_coeff_free(my_coeffs(ii))
1095  end do
1096 
1097  ABI_DATATYPE_DEALLOCATE(my_coeffs)
1098  ABI_DEALLOCATE(buffsize)
1099  ABI_DEALLOCATE(buffdisp)
1100  ABI_DEALLOCATE(buffGF)
1101  ABI_DEALLOCATE(coeff_values)
1102  ABI_DEALLOCATE(energy_coeffs)
1103  ABI_DEALLOCATE(energy_coeffs_tmp)
1104  ABI_DEALLOCATE(fcart_coeffs)
1105  ABI_DEALLOCATE(fcart_coeffs_tmp)
1106  ABI_DEALLOCATE(gf_mpi)
1107  ABI_DEALLOCATE(gf_values)
1108  ABI_DEALLOCATE(list_coeffs)
1109  ABI_DEALLOCATE(list_coeffs_tmp)
1110  ABI_DEALLOCATE(list_coeffs_tmp2)
1111  ABI_DEALLOCATE(my_coeffindexes)
1112  ABI_DEALLOCATE(my_coefflist)
1113  ABI_DEALLOCATE(singular_coeffs)
1114  ABI_DEALLOCATE(strten_coeffs)
1115  ABI_DEALLOCATE(strten_coeffs_tmp)
1116  ABI_DEALLOCATE(stat_coeff)
1117 
1118 end subroutine fit_polynomial_coeff_fit

m_fit_polynomial_coeff/fit_polynomial_coeff_getCoeffBound [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_getCoeffBound

FUNCTION

 This routine fit a list of possible model.
 Return in the isPositive array:

INPUTS

 NEED TO UPDATE
 eff_pot<type(effective_potential)> = effective potential
 hist<type(abihist)> = The history of the MD (or snapshot of DFT
 comm = MPI communicator
 verbose  = optional, flag for the verbose mode

OUTPUT

PARENTS

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

1361 subroutine fit_polynomial_coeff_getCoeffBound(eff_pot,coeffs_out,hist,ncoeff_bound,comm,verbose)
1362 
1363 
1364 !This section has been created automatically by the script Abilint (TD).
1365 !Do not modify the following lines by hand.
1366 #undef ABI_FUNC
1367 #define ABI_FUNC 'fit_polynomial_coeff_getCoeffBound'
1368 !End of the abilint section
1369 
1370  implicit none
1371 
1372 !Arguments ------------------------------------
1373  !scalars
1374  integer,intent(in) :: comm
1375  integer,intent(out) :: ncoeff_bound
1376  logical,optional,intent(in) :: verbose
1377 !arrays
1378  type(abihist),intent(inout) :: hist
1379  type(effective_potential_type),target,intent(inout) :: eff_pot
1380  type(polynomial_coeff_type),allocatable,intent(out) :: coeffs_out(:)
1381 !Local variables-------------------------------
1382 !scalar
1383  integer :: counter,icoeff,icoeff_bound,idisp,istrain,ii
1384  integer :: istart,iterm,ndisp,nstrain,nterm,ncoeff_model,ncoeff_in,ncoeff_max
1385  real(dp):: weight
1386  logical :: need_verbose
1387 !arrays
1388  integer,allocatable :: atindx(:,:),cells(:,:,:),direction(:)
1389  integer,allocatable :: power_disps(:),power_strain(:),strain(:)
1390  type(polynomial_term_type),dimension(:),allocatable :: terms
1391  integer,allocatable :: odd_coeff(:),need_bound(:)
1392  type(polynomial_coeff_type),pointer :: coeffs_in(:)
1393  type(polynomial_coeff_type),allocatable :: coeffs_test(:)
1394  character(len=5),allocatable :: symbols(:)
1395  character(len=200):: name
1396  character(len=500) :: msg
1397 ! *************************************************************************
1398 
1399 
1400 !set the inputs varaibles
1401  ncoeff_model =  eff_pot%anharmonics_terms%ncoeff
1402  coeffs_in => eff_pot%anharmonics_terms%coefficients
1403 
1404 !Do check
1405  if(ncoeff_model == 0)then
1406    write(msg,'(a)')'ncoeff_model must be different to 0'
1407    MSG_BUG(msg)
1408  end if
1409 
1410 !Map the hist in order to be consistent with the supercell into reference_effective_potential
1411  call effective_potential_file_mapHistToRef(eff_pot,hist,comm)
1412 
1413 !Initialisation of optional arguments
1414  need_verbose = .TRUE.
1415  if(present(verbose)) need_verbose = verbose
1416 
1417  write(msg, '(a)' ) ' Detection of the unbound coefficients'
1418  if(need_verbose)call wrtout(std_out,msg,'COLL')
1419 
1420 !Allocation
1421  ncoeff_max = 2 * ncoeff_model
1422  ABI_ALLOCATE(odd_coeff,(ncoeff_max))
1423  ABI_ALLOCATE(need_bound,(ncoeff_max))
1424 
1425  ABI_ALLOCATE(symbols,(eff_pot%crystal%natom))
1426  call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,&
1427 &                     symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
1428 
1429 
1430  ABI_DATATYPE_ALLOCATE(coeffs_test,(ncoeff_max))
1431 
1432  do icoeff=1,ncoeff_model
1433    call polynomial_coeff_init(coeffs_in(icoeff)%coefficient,coeffs_in(icoeff)%nterm,&
1434 &                             coeffs_test(icoeff),coeffs_in(icoeff)%terms,&
1435 &                             coeffs_in(icoeff)%name,check=.false.)
1436  end do
1437 
1438 !array to know which coeff has to be bound
1439  need_bound(:) = 1
1440  counter = 0
1441  ncoeff_in = ncoeff_model
1442 
1443  do while(.not.all(need_bound == 0).and.counter<1)
1444 !  Get the coefficients with odd coefficient
1445    odd_coeff = 0
1446    if(counter>0) then
1447      need_bound(1:ncoeff_in) = 0
1448      icoeff_bound = ncoeff_in
1449    else
1450      icoeff_bound = 1
1451    end if
1452 
1453    do icoeff=icoeff_bound,ncoeff_model
1454      if(any(mod(coeffs_in(icoeff)%terms(1)%power_disp(:),2)/=0))then
1455        odd_coeff(icoeff) = 1
1456      end if
1457      if(any(mod(coeffs_in(icoeff)%terms(1)%power_strain(:),2)/=0))then
1458        odd_coeff(icoeff) = 1
1459      end if
1460      if(odd_coeff(icoeff) == 0 .and. coeffs_in(icoeff)%coefficient > zero) then
1461        need_bound(icoeff) = 0
1462      else
1463         need_bound(icoeff) = 1
1464      end if
1465    end do
1466    if(need_verbose)then
1467      write(msg, '(a)' ) ' The following coefficients need to be bound:'
1468      call wrtout(std_out,msg,'COLL')
1469      do icoeff=1,ncoeff_model
1470        if(need_bound(icoeff) == 1)then
1471          write(msg, '(2a)' ) ' =>',trim(coeffs_in(icoeff)%name)
1472          call wrtout(std_out,msg,'COLL')
1473        end if
1474      end do
1475    end if
1476 
1477 
1478    icoeff_bound = ncoeff_in + 1
1479    if(counter==0)then
1480      istart = 1
1481    else
1482      istart = ncoeff_in
1483    end if
1484 
1485    ncoeff_bound = count(need_bound(istart:ncoeff_model)==1)
1486 
1487    do icoeff=istart,ncoeff_model
1488      if(need_bound(icoeff)==1)then
1489 
1490        nterm = coeffs_in(icoeff)%nterm
1491        ndisp = coeffs_in(icoeff)%terms(1)%ndisp
1492        nstrain = coeffs_in(icoeff)%terms(1)%nstrain
1493 
1494        ABI_ALLOCATE(terms,(nterm))
1495        ABI_ALLOCATE(atindx,(2,ndisp))
1496        ABI_ALLOCATE(cells,(3,2,ndisp))
1497        ABI_ALLOCATE(direction,(ndisp))
1498        ABI_ALLOCATE(power_disps,(ndisp))
1499        ABI_ALLOCATE(power_strain,(nstrain))
1500        ABI_ALLOCATE(strain,(nstrain))
1501 
1502        do iterm=1,coeffs_in(icoeff)%nterm
1503          atindx(:,:) = coeffs_in(icoeff)%terms(iterm)%atindx(:,:)
1504          cells(:,:,:) = coeffs_in(icoeff)%terms(iterm)%cell(:,:,:)
1505          direction(:) = coeffs_in(icoeff)%terms(iterm)%direction(:)
1506          power_strain(:) = coeffs_in(icoeff)%terms(iterm)%power_strain(:)
1507          power_disps(:) = coeffs_in(icoeff)%terms(iterm)%power_disp(:)
1508          strain(:) =  coeffs_in(icoeff)%terms(iterm)%strain(:)
1509          weight =  1
1510          do idisp=1,ndisp
1511            if(mod(power_disps(idisp),2) /= 0) then
1512              power_disps(idisp) = power_disps(idisp) + 1
1513            else
1514              power_disps(idisp) = power_disps(idisp) + 2
1515            end if
1516          end do
1517          do istrain=1,nstrain
1518            if(mod(power_strain(istrain),2) /= 0)then
1519              power_strain(istrain) = power_strain(istrain) + 1
1520            else
1521              if(power_strain(istrain) < 4 ) power_strain(istrain) = power_strain(istrain) + 2
1522            end if
1523          end do
1524 
1525          call polynomial_term_init(atindx,cells,direction,ndisp,nstrain,terms(iterm),&
1526 &                                  power_disps,power_strain,strain,weight,check=.true.)
1527        end do
1528 
1529        name = ""
1530        call polynomial_coeff_init(one,nterm,coeffs_test(icoeff_bound),terms,name,check=.true.)
1531        call polynomial_coeff_getName(name,coeffs_test(icoeff_bound),symbols,recompute=.TRUE.)
1532        call polynomial_coeff_SetName(name,coeffs_test(icoeff_bound))
1533 
1534 !      Deallocate the terms
1535        do iterm=1,nterm
1536          call polynomial_term_free(terms(iterm))
1537        end do
1538        ABI_DEALLOCATE(terms)
1539        ABI_DEALLOCATE(atindx)
1540        ABI_DEALLOCATE(cells)
1541        ABI_DEALLOCATE(direction)
1542        ABI_DEALLOCATE(power_disps)
1543        ABI_DEALLOCATE(power_strain)
1544        ABI_DEALLOCATE(strain)
1545 
1546        icoeff_bound = icoeff_bound  + 1
1547 
1548      end if
1549    end do
1550 
1551 
1552    if(counter==0)ncoeff_model = ncoeff_model + ncoeff_bound
1553 !   call effective_potential_setCoeffs(coeffs_test,eff_pot,ncoeff_model)
1554    call fit_polynomial_coeff_fit(eff_pot,(/0/),(/0/),hist,0,(/0,0/),1,0,&
1555 &             -1,1,comm,verbose=.true.,positive=.false.)
1556 
1557    coeffs_in => eff_pot%anharmonics_terms%coefficients
1558 
1559    counter = counter + 1
1560  end do
1561 
1562  ABI_ALLOCATE(coeffs_out,(ncoeff_bound))
1563  do ii=1,ncoeff_bound
1564    icoeff_bound = ncoeff_in + ii
1565    call polynomial_coeff_init(one,coeffs_test(icoeff_bound)%nterm,coeffs_out(ii),&
1566 &                    coeffs_test(icoeff_bound)%terms,coeffs_test(icoeff_bound)%name,check=.true.)
1567  end do
1568 !Deallocation
1569  do ii=ncoeff_model,ncoeff_max
1570    call polynomial_coeff_free(coeffs_test(ii))
1571  end do
1572 
1573 !Deallocation
1574  do icoeff=1,ncoeff_max
1575    call polynomial_coeff_free(coeffs_test(icoeff))
1576  end do
1577  ABI_DATATYPE_DEALLOCATE(coeffs_test)
1578  ABI_DEALLOCATE(odd_coeff)
1579  ABI_DEALLOCATE(need_bound)
1580  ABI_DEALLOCATE(symbols)
1581 
1582 
1583 end subroutine fit_polynomial_coeff_getCoeffBound

m_fit_polynomial_coeff/fit_polynomial_coeff_getFS [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_getFS

FUNCTION

 Compute all the matrix elements of eq.11 and 12 in PRB95,094115 (2017) [[cite:Escorihuela-Sayalero2017]]

INPUTS

 coefficients(ncoeff)          = type(polynomial_coeff_type)
 du_delta(6,3,natom_sc,ntime)  = Variation to displacements wrt to the strain (Bohr)
 displacement(3,natom_sc,ntime)= Atomic displacement wrt to the reference (Bohr)
 natom_sc = Number of atoms in the supercell
 natom_uc = Number of atoms in the unit cell
 ncoeff = Number of coefficients
 ntime = Number of time in the history
 sc_size(3) = Size of the supercell
 strain(6,ntime) = Strain
 ucvol(ntime) = Volume of the supercell for each time (Bohr^3)
 cells(ncell) = Indexes of the cell treat by this CPU
 ncell = Number of cell treat by this CPU
 index_cells(ncell,3) = Indexes of the cells (1 1 1, 0 0 0 for instance) treat by this CPU
 comm  = MPI communicator

OUTPUT

 fcart_out(ncoeff,3,natom,ntime) = value of the forces for each coefficient
                                   (-1 factor is taking into acount) (Ha/Bohr)
 strten_out(ncoeff,3,natom,ntime)= value of the stresses for each coefficient
                                   (-1/ucvol factor is taking into acount) (Ha/Bohr^3)
 energy_out(ncoeff,ntime)        = value of the energy for each  coefficient (Ha)

PARENTS

      m_fit_polynomial_coeff

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

1965 subroutine fit_polynomial_coeff_getFS(coefficients,du_delta,displacement,energy_out,fcart_out,&
1966 &                                     natom_sc,natom_uc,ncoeff_max,ntime,sc_size,strain,strten_out,&
1967 &                                     ucvol,coeffs,ncoeff)
1968 
1969 
1970 !This section has been created automatically by the script Abilint (TD).
1971 !Do not modify the following lines by hand.
1972 #undef ABI_FUNC
1973 #define ABI_FUNC 'fit_polynomial_coeff_getFS'
1974 !End of the abilint section
1975 
1976  implicit none
1977 
1978 !Arguments ------------------------------------
1979 !scalars
1980  integer,intent(in) :: natom_sc,natom_uc,ncoeff_max,ntime
1981  integer,intent(in) :: ncoeff
1982 !arrays
1983  integer,intent(in) :: sc_size(3)
1984  integer,intent(in) :: coeffs(ncoeff_max)
1985  real(dp),intent(in) :: du_delta(6,3,natom_sc,ntime)
1986  real(dp),intent(in) :: displacement(3,natom_sc,ntime)
1987  real(dp),intent(in) :: strain(6,ntime),ucvol(ntime)
1988  real(dp),intent(out):: energy_out(ncoeff,ntime)
1989  real(dp),intent(out) :: fcart_out(3,natom_sc,ncoeff,ntime)
1990  real(dp),intent(out) :: strten_out(6,ntime,ncoeff)
1991  type(polynomial_coeff_type), intent(in) :: coefficients(ncoeff_max)
1992 !Local variables-------------------------------
1993 !scalar
1994  integer :: i1,i2,i3,ia1,ia2,ib1,ib2,ii,icell,icoeff,icoeff_tmp
1995  integer :: idir1,idir2,idisp1,idisp2,idisp1_strain,idisp2_strain
1996  integer :: iterm,itime,ndisp,ndisp_tot,nstrain,power_disp,power_strain
1997  real(dp):: disp1,disp2,tmp1,tmp2,tmp3,weight
1998 !arrays
1999  integer :: cell_atoma1(3),cell_atoma2(3)
2000  integer :: cell_atomb1(3),cell_atomb2(3)
2001 
2002 ! *************************************************************************
2003 
2004 
2005 !1-Get forces and stresses from the model
2006 !  Initialisation of variables
2007  fcart_out(:,:,:,:) = zero
2008  strten_out(:,:,:)  = zero
2009  energy_out(:,:)    = zero
2010 
2011  icell = 0; ib1=0; ia1=0
2012  do i1=1,sc_size(1)
2013    do i2=1,sc_size(2)
2014      do i3=1,sc_size(3)
2015        ii = icell*natom_uc
2016        icell = icell + 1
2017 !      Loop over configurations
2018        do itime=1,ntime
2019 !       Loop over coefficients
2020          do icoeff_tmp=1,ncoeff
2021            icoeff = coeffs(icoeff_tmp)
2022 !          Loop over terms of this coefficient
2023            do iterm=1,coefficients(icoeff)%nterm
2024              ndisp = coefficients(icoeff)%terms(iterm)%ndisp
2025              nstrain = coefficients(icoeff)%terms(iterm)%nstrain
2026              ndisp_tot = ndisp + nstrain
2027 !            Set the weight of this term
2028              weight =coefficients(icoeff)%terms(iterm)%weight
2029              tmp1 = one
2030 !            Loop over displacement and strain
2031              do idisp1=1,ndisp_tot
2032 
2033 !              Set to one the acculation of forces and strain
2034                tmp2 = one
2035                tmp3 = one
2036 !              Strain case idir => -6, -5, -4, -3, -2 or -1
2037                if (idisp1 > ndisp)then
2038                  idisp1_strain = idisp1 - ndisp
2039                  power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain)
2040 !                Get the direction of the displacement or strain
2041                  idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain)
2042                  if(abs(strain(idir1,itime)) > tol10)then
2043 !                  Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
2044                    tmp1 = tmp1 * (strain(idir1,itime))**power_strain
2045                    if(power_strain > 1) then
2046 !                    Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
2047                      tmp3 = tmp3 *  power_strain*(strain(idir1,itime))**(power_strain-1)
2048                    end if
2049                  else
2050                    tmp1 = zero
2051                    if(power_strain > 1) then
2052                      tmp3 = zero
2053                    end if
2054                  end if
2055                else
2056 !                Set the power_disp of the displacement:
2057                  power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1)
2058 !                Get the direction of the displacement or strain
2059                  idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1)
2060 !                Displacement case idir = 1, 2  or 3
2061 !                indexes of the cell of the atom a
2062                  cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1)
2063                  if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then
2064 !                  if the cell is not 0 0 0 we apply PBC:
2065                    cell_atoma1(1) =  i1 + cell_atoma1(1)
2066                    cell_atoma1(2) =  i2 + cell_atoma1(2)
2067                    cell_atoma1(3) =  i3 + cell_atoma1(3)
2068                    call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3))
2069 !                  index of the first atom (position in the supercell if the cell is not 0 0 0)
2070                    ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
2071 &                        (cell_atoma1(2)-1)*sc_size(3)*natom_uc+&
2072 &                        (cell_atoma1(3)-1)*natom_uc+&
2073 &                        coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
2074                  else
2075 !                  index of the first atom (position in the supercell if the cell is 0 0 0)
2076                    ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
2077                  end if
2078 
2079 !                indexes of the cell of the atom b  (with PBC) same as ia1
2080                  cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1)
2081                  if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then
2082                    cell_atomb1(1) =  i1 + cell_atomb1(1)
2083                    cell_atomb1(2) =  i2 + cell_atomb1(2)
2084                    cell_atomb1(3) =  i3 + cell_atomb1(3)
2085                    call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3))
2086 
2087 !                  index of the second atom in the (position in the supercell  if the cell is not 0 0 0)
2088                    ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
2089 &                        (cell_atomb1(2)-1)*sc_size(3)*natom_uc+&
2090 &                        (cell_atomb1(3)-1)*natom_uc+&
2091 &                        coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
2092                  else
2093 !                  index of the first atom (position in the supercell if the cell is 0 0 0)
2094                    ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
2095                  end if
2096 
2097 !                Get the displacement for the both atoms
2098                  disp1 = displacement(idir1,ia1,itime)
2099                  disp2 = displacement(idir1,ib1,itime)
2100 
2101                  if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then
2102 !                  Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
2103                    tmp1 = tmp1 * (disp1-disp2)**power_disp
2104                    if(power_disp > 1) then
2105 !                    Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...))
2106                      tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1)
2107                    end if
2108                  else
2109                    tmp1 = zero
2110                    if(power_disp > 1) then
2111                      tmp2 = zero
2112                    end if
2113                  end if
2114                end if
2115 
2116                do idisp2=1,ndisp_tot
2117                  if(idisp2 /= idisp1) then
2118 
2119 !                  Strain case
2120                    if (idisp2 > ndisp)then
2121                      idisp2_strain = idisp2 - ndisp
2122                      idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain)
2123 !                    Set the power_strain of the strain:
2124                      power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain)
2125 !                    Accumulate energy forces
2126                      tmp2 = tmp2 * (strain(idir2,itime))**power_strain
2127 !                    Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
2128                      tmp3 = tmp3 * (strain(idir2,itime))**power_strain
2129 !                  Atomic displacement case
2130                    else
2131 !                    Set the power_disp of the displacement:
2132                      power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2)
2133 !                    Set the direction of the displacement:
2134                      idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2)
2135 
2136                      cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2)
2137                      if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then
2138                        cell_atoma2(1) =  i1 + cell_atoma2(1)
2139                        cell_atoma2(2) =  i2 + cell_atoma2(2)
2140                        cell_atoma2(3) =  i3 + cell_atoma2(3)
2141                        call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3))
2142 !                      index of the first atom (position in the supercell and direction)
2143 !                      if the cell of the atom a is not 0 0 0 (may happen)
2144                        ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
2145 &                            (cell_atoma2(2)-1)*sc_size(3)*natom_uc+&
2146 &                            (cell_atoma2(3)-1)*natom_uc+&
2147 &                        coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
2148                      else
2149 !                      index of the first atom (position in the supercell and direction)
2150                        ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
2151                      end if
2152 
2153                      cell_atomb2 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2)
2154 
2155                      if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then
2156 !                      indexes of the cell2 (with PBC)
2157                        cell_atomb2(1) =  i1 + cell_atomb2(1)
2158                        cell_atomb2(2) =  i2 + cell_atomb2(2)
2159                        cell_atomb2(3) =  i3 + cell_atomb2(3)
2160                        call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3))
2161 
2162 !                      index of the second atom in the (position in the supercell)
2163                        ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
2164 &                            (cell_atomb2(2)-1)*sc_size(3)*natom_uc+&
2165 &                            (cell_atomb2(3)-1)*natom_uc+&
2166 &                            coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
2167                      else
2168                        ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
2169                      end if
2170 
2171                      disp1 = displacement(idir2,ia2,itime)
2172                      disp2 = displacement(idir2,ib2,itime)
2173 
2174                      tmp2 = tmp2 * (disp1-disp2)**power_disp
2175                      tmp3 = tmp3 * (disp1-disp2)**power_disp
2176 
2177                    end if
2178                  end if
2179                end do
2180 
2181                if(idisp1 > ndisp)then
2182 !                Accumule stress tensor
2183                  strten_out(idir1,itime,icoeff_tmp) = strten_out(idir1,itime,icoeff_tmp) + &
2184 &                                                      weight * tmp3 / ucvol(itime)
2185                else
2186 !                Accumule  forces
2187                  fcart_out(idir1,ia1,icoeff_tmp,itime)=fcart_out(idir1,ia1,icoeff_tmp,itime)+weight*tmp2
2188                  fcart_out(idir1,ib1,icoeff_tmp,itime)=fcart_out(idir1,ib1,icoeff_tmp,itime)-weight*tmp2
2189                end if
2190              end do
2191 
2192 !            accumule energy
2193              energy_out(icoeff_tmp,itime) = energy_out(icoeff_tmp,itime) +  weight * tmp1
2194 
2195            end do!End do iterm
2196          end do!End do coeff
2197        end do!End time
2198      end do!End do i3
2199    end do!End do i2
2200  end do!End do i1
2201 
2202 !ADD variation of the atomic displacement due to the strain
2203  do icoeff=1,ncoeff
2204    do itime=1,ntime
2205      do ia1=1,natom_sc
2206        do idir1=1,3
2207          do idir2=1,6
2208            strten_out(idir2,itime,icoeff) = strten_out(idir2,itime,icoeff) + &
2209 &                     du_delta(idir2,idir1,ia1,itime)*fcart_out(idir1,ia1,icoeff,itime)/ucvol(itime)
2210          end do
2211        end do
2212      end do
2213    end do
2214  end do
2215 
2216 ! multiply by -1
2217  fcart_out(:,:,:,:) = -1 * fcart_out(:,:,:,:)
2218 
2219 end subroutine fit_polynomial_coeff_getFS

m_fit_polynomial_coeff/fit_polynomial_coeff_getPositive [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_getPositive

FUNCTION

 This routine fit a list of possible model.
 Return in the isPositive array:
   0 if the model ii does not contain possive coefficients
   1 if the model ii contain possive coefficients

INPUTS

 eff_pot<type(effective_potential)> = effective potential
 hist<type(abihist)> = The history of the MD (or snapshot of DFT
 coeff_values(nmodel,ncoeff) = values of the coefficients for each model
 isPositive(nmodel) = see description below
 list_coeff(nmodel,ncoeff) = list of the models
 ncoeff = number of coeff per model
 nfixcoeff = will not test the nfixcoeff first coeffcients
 nmodel = number of model
 comm = MPI communicator
 verbose  = optional, flag for the verbose mode

OUTPUT

 eff_pot = effective potential datatype with new fitted coefficients

PARENTS

      mover_effpot

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

1155 subroutine fit_polynomial_coeff_getPositive(eff_pot,hist,coeff_values,isPositive,list_coeff,ncoeff,&
1156 &                                           nfixcoeff,nmodel,comm,verbose)
1157 
1158 
1159 !This section has been created automatically by the script Abilint (TD).
1160 !Do not modify the following lines by hand.
1161 #undef ABI_FUNC
1162 #define ABI_FUNC 'fit_polynomial_coeff_getPositive'
1163 !End of the abilint section
1164 
1165  implicit none
1166 
1167 !Arguments ------------------------------------
1168 !scalars
1169  integer,intent(in) :: ncoeff,nfixcoeff,nmodel,comm
1170 !arrays
1171  integer,intent(in)  :: list_coeff(nmodel,ncoeff)
1172  integer,intent(out) :: isPositive(nmodel)
1173  real(dp),intent(out) :: coeff_values(nmodel,ncoeff)
1174  type(effective_potential_type),intent(inout) :: eff_pot
1175  type(abihist),intent(inout) :: hist
1176  logical,optional,intent(in) :: verbose
1177 !Local variables-------------------------------
1178 !scalar
1179  integer :: ierr,ii,info,imodel,my_nmodel,nmodel_alone
1180  integer :: master,my_rank,ncoeff_tot,natom_sc,ncell
1181  integer :: nproc,ntime
1182  logical :: iam_master,need_verbose
1183 !arrays
1184  integer :: sc_size(3)
1185  integer,allocatable  :: list_coeffs(:),my_modelindexes(:),my_modellist(:)
1186  real(dp),allocatable :: energy_coeffs(:,:),fcart_coeffs(:,:,:,:), strten_coeffs(:,:,:)
1187  type(polynomial_coeff_type),allocatable :: coeffs_in(:)
1188  type(fit_data_type) :: fit_data
1189  character(len=500) :: message
1190 ! *************************************************************************
1191 
1192 !MPI variables
1193  master = 0
1194  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1195  iam_master = (my_rank == master)
1196 
1197 !Initialisation of optional arguments
1198  need_verbose = .TRUE.
1199  if(present(verbose)) need_verbose = verbose
1200 
1201 !Get the list of coefficients from the eff_pot
1202  if(eff_pot%anharmonics_terms%ncoeff > 0)then
1203 !  Copy the initial coefficients array
1204    ncoeff_tot = eff_pot%anharmonics_terms%ncoeff
1205    ABI_DATATYPE_ALLOCATE(coeffs_in,(ncoeff_tot))
1206    do ii=1,ncoeff_tot
1207      call polynomial_coeff_init(eff_pot%anharmonics_terms%coefficients(ii)%coefficient,&
1208 &                               eff_pot%anharmonics_terms%coefficients(ii)%nterm,&
1209 &                               coeffs_in(ii),&
1210 &                               eff_pot%anharmonics_terms%coefficients(ii)%terms,&
1211 &                               eff_pot%anharmonics_terms%coefficients(ii)%name,&
1212 &                               check=.false.)
1213    end do
1214  end if
1215 
1216 !Reset the output (we free the memory)
1217  call effective_potential_freeCoeffs(eff_pot)
1218 
1219 !if the number of atoms in reference supercell into effpot is not corret,
1220 !wrt to the number of atom in the hist, we set map the hist and set the good
1221 !supercell
1222  if (size(hist%xred,2) /= eff_pot%supercell%natom) then
1223    call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose)
1224  end if
1225 
1226 !Initialisation of constants
1227  natom_sc   = eff_pot%supercell%natom
1228  ncell      = eff_pot%supercell%ncells
1229  ntime      = hist%mxhist
1230  do ii = 1, 3
1231    sc_size(ii) = eff_pot%supercell%rlatt(ii,ii)
1232  end do
1233 
1234 !Initialisation of arrays:
1235  ABI_ALLOCATE(list_coeffs,(ncoeff_tot))
1236  list_coeffs  = 0
1237  do ii = 1,ncoeff_tot
1238    list_coeffs(ii) = ii
1239  end do
1240 
1241 !Get the decomposition for each coefficients of the forces and stresses for
1242 !each atoms and each step  equations 11 & 12 of  PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
1243  if(need_verbose)then
1244    write(message, '(a)' ) ' Initialisation of the fit process...'
1245    call wrtout(std_out,message,'COLL')
1246  end if
1247 !Before the fit, compute constants with fit_data_compute.
1248 !Conpute the strain of each configuration.
1249 !Compute the displacmeent of each configuration.
1250 !Compute the variation of the displacement due to strain of each configuration.
1251 !Compute fixed forces and stresse and get the standard deviation.
1252 !Compute Sheppard and al Factors  \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
1253  call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose)
1254 
1255 !Get the decomposition for each coefficients of the forces,stresses and energy for
1256 !each atoms and each step  (see equations 11 & 12 of  
1257 ! PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]] + allocation
1258  ABI_ALLOCATE(energy_coeffs,(ncoeff_tot,ntime))
1259  ABI_ALLOCATE(fcart_coeffs,(3,natom_sc,ncoeff_tot,ntime))
1260  ABI_ALLOCATE(strten_coeffs,(6,ntime,ncoeff_tot))
1261 
1262  call fit_polynomial_coeff_getFS(coeffs_in,fit_data%training_set%du_delta,&
1263 &                                fit_data%training_set%displacement,&
1264 &                                energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
1265 &                                ncoeff_tot,ntime,sc_size,&
1266 &                                fit_data%training_set%strain,strten_coeffs,&
1267 &                                fit_data%training_set%ucvol,list_coeffs,ncoeff_tot)
1268 
1269 
1270 !set MPI, really basic stuff...
1271  nmodel_alone = mod(nmodel,nproc)
1272  my_nmodel = int(aint(real(nmodel,sp)/(nproc)))
1273 
1274  if(my_rank >= (nproc-nmodel_alone)) then
1275    my_nmodel = my_nmodel  + 1
1276  end if
1277 
1278  ABI_ALLOCATE(my_modelindexes,(my_nmodel))
1279  ABI_ALLOCATE(my_modellist,(my_nmodel))
1280 
1281 !2:compute the number of model and the list of the corresponding for each CPU.
1282  do imodel=1,my_nmodel
1283    if(my_rank >= (nproc-nmodel_alone))then
1284      my_modelindexes(imodel)=(int(aint(real(nmodel,sp)/nproc)))*(my_rank)+&
1285 &                              (my_rank - (nproc-nmodel_alone)) + imodel
1286      my_modellist(imodel) = imodel
1287    else
1288      my_modelindexes(imodel)=(my_nmodel)*(my_rank)  + imodel
1289      my_modellist(imodel) = imodel
1290   end if
1291  end do
1292 
1293 
1294 !Start fit process
1295  isPositive   = 0
1296  coeff_values = zero
1297  do ii=1,my_nmodel
1298    imodel = my_modelindexes(ii)
1299    call fit_polynomial_coeff_solve(coeff_values(imodel,1:ncoeff),fcart_coeffs,fit_data%fcart_diff,&
1300 &                                  energy_coeffs,fit_data%energy_diff,info,&
1301 &                                  list_coeff(imodel,1:ncoeff),natom_sc,ncoeff,&
1302 &                                  ncoeff_tot,ntime,strten_coeffs,fit_data%strten_diff,&
1303 &                                  fit_data%training_set%sqomega)
1304 
1305    if(info==0)then
1306 
1307      if (any(coeff_values(imodel,nfixcoeff+1:ncoeff) < zero))then
1308 !       coeff_values(imodel,:) = zero
1309        isPositive(imodel) = 0
1310      else
1311        isPositive(imodel) = 1
1312      end if
1313    end if
1314  end do
1315 
1316  call xmpi_sum(isPositive, comm, ierr)
1317  call xmpi_sum(coeff_values, comm, ierr)
1318 
1319 !Deallocation of arrays
1320  do ii=1,ncoeff_tot
1321    call polynomial_coeff_free(coeffs_in(ii))
1322  end do
1323  call fit_data_free(fit_data)
1324  ABI_DATATYPE_DEALLOCATE(coeffs_in)
1325  ABI_DEALLOCATE(energy_coeffs)
1326  ABI_DEALLOCATE(fcart_coeffs)
1327  ABI_DEALLOCATE(list_coeffs)
1328  ABI_DEALLOCATE(my_modelindexes)
1329  ABI_DEALLOCATE(my_modellist)
1330  ABI_DEALLOCATE(strten_coeffs)
1331 
1332 end subroutine fit_polynomial_coeff_getPositive

m_fit_polynomial_coeff/fit_polynomial_coeff_solve [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_coeff_solve

FUNCTION

 Build and the solve the system to get the values of the coefficients
 This routine solves the linear system proposed by 
 C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]

INPUTS

 fcart_coeffs(3,natom_sc,ncoeff_max,ntime) = List of the values of the contribution to the
                                             cartesian forces for all coefficients
                                             for each direction and each time
 fcart_diff(3,natom,ntime) = Difference of cartesian forces between DFT calculation and
                             fixed part of the model (more often harmonic part)
 energy_coeffs(ncoeff,ntime)   = value of the energy for each  coefficient (Ha)
 energy_diff(ntime) = Difference of energ ybetween DFT calculation and fixed part
                             of the model (more often harmonic part)
 list_coeffs(ncoeff_fit) = List with the index of the coefficients used for this model
 natom = Number of atoms
 ncoeff_fit = Number of coeff for the fit (dimension of the system)
 ncoeff_max = Maximum number of coeff in the list
 ntime = Number of time (number of snapshot, number of md step...)
 strten_coeffs(6,ntime,ncoeff_max) = List of the values of the contribution to the stress tensor
                                      of  the coefficients for each direction,time
 strten_diff(6,natom) = Difference of stress tensor between DFT calculation and
                        fixed part of the model (more often harmonic part)
 sqomega(ntime) =  Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]

OUTPUT

 coefficients(ncoeff_fit) = Values of the coefficients
 info_out = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, U(i,i) computed in DOUBLE PRECISION is
                exactly zero.  The factorization has been completed,
                but the factor U is exactly singular, so the solution
                could not be computed.  = 0:  successful exit
          information from the subroutine dsgesv in LAPACK

PARENTS

      m_fit_polynomial_coeff

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

1635 subroutine fit_polynomial_coeff_solve(coefficients,fcart_coeffs,fcart_diff,energy_coeffs,energy_diff,&
1636 &                                     info_out,list_coeffs,natom,ncoeff_fit,ncoeff_max,ntime,&
1637 &                                     strten_coeffs,strten_diff,sqomega)
1638 
1639 
1640 !This section has been created automatically by the script Abilint (TD).
1641 !Do not modify the following lines by hand.
1642 #undef ABI_FUNC
1643 #define ABI_FUNC 'fit_polynomial_coeff_solve'
1644 !End of the abilint section
1645 
1646  implicit none
1647 
1648 !Arguments ------------------------------------
1649 !scalars
1650  integer,intent(in)  :: natom,ncoeff_fit,ncoeff_max,ntime
1651  integer,intent(out) :: info_out
1652  !arrays
1653  real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime)
1654  real(dp),intent(in) :: energy_diff(ntime)
1655  integer,intent(in)  :: list_coeffs(ncoeff_fit)
1656  real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime)
1657  real(dp),intent(in) :: fcart_diff(3,natom,ntime)
1658  real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max)
1659  real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime)
1660  real(dp),intent(out):: coefficients(ncoeff_fit)
1661 !Local variables-------------------------------
1662 !scalar
1663  integer :: ia,itime,icoeff,jcoeff,icoeff_tmp,jcoeff_tmp,mu,LDA,LDB,LDX,LDAF,N,NRHS
1664  real(dp):: efact,ffact,sfact,ftmpA,stmpA,ftmpB,stmpB,etmpA,etmpB,fmu,fnu,smu,snu,emu,enu
1665  integer :: INFO,ITER
1666  real(dp):: RCOND
1667  real(dp):: fcart_coeffs_tmp(3,natom,ntime)
1668  real(dp),allocatable:: AF(:,:),BERR(:),FERR(:),WORK(:),C(:),R(:)
1669  integer,allocatable :: IPIV(:),IWORK(:),SWORK(:)
1670 !arrays
1671  real(dp),allocatable :: A(:,:),B(:,:)
1672  character(len=1) :: FACT,EQUED,TRANS
1673 ! character(len=500) :: message
1674 ! *************************************************************************
1675 
1676 !0-Set variables for the
1677  N    = ncoeff_fit; NRHS = 1; LDA  = ncoeff_fit; LDB  = ncoeff_fit; LDX  = ncoeff_fit
1678  LDAF = ncoeff_fit;  RCOND = zero; INFO  = 0; TRANS='N'; EQUED='N'; FACT='N'
1679 
1680 !Set the factors
1681  ffact = one/(3*natom*ntime)
1682  sfact = one/(6*ntime)
1683  efact = one/(ntime)
1684 
1685 !0-Allocation
1686  ABI_ALLOCATE(A,(LDA,N))
1687  ABI_ALLOCATE(B,(LDB,NRHS))
1688  ABI_ALLOCATE(AF,(LDAF,N))
1689  ABI_ALLOCATE(IPIV,(N))
1690  ABI_ALLOCATE(R,(N))
1691  ABI_ALLOCATE(C,(N))
1692  ABI_ALLOCATE(FERR,(NRHS))
1693  ABI_ALLOCATE(BERR,(NRHS))
1694  ABI_ALLOCATE(WORK,(4*N))
1695  ABI_ALLOCATE(IWORK,(N))
1696  ABI_ALLOCATE(SWORK,(N*(N+NRHS)))
1697  A=zero; B=zero;
1698  AF = zero; IPIV = 1;
1699  R = one; C = one;
1700  FERR = zero; BERR = zero
1701  IWORK = 0; WORK = 0
1702 
1703 !1-Get forces and stresses from the model and fill A
1704 !  Fill alsor B with the forces and stresses from
1705 !  the DFT snapshot and the model
1706 !  See equation 17 of PRB95 094115 (2017) [[cite:Escorihuela-Sayalero2017]]
1707  do icoeff=1,ncoeff_fit
1708    icoeff_tmp = list_coeffs(icoeff)
1709    fcart_coeffs_tmp(:,:,:) = fcart_coeffs(:,:,icoeff_tmp,:)
1710    ftmpA= zero; ftmpB = zero
1711    stmpA= zero; stmpB = zero
1712    etmpA= zero; etmpB = zero
1713 !  loop over the configuration
1714    do itime=1,ntime
1715 !    Fill energy
1716      emu = energy_coeffs(icoeff_tmp,itime)
1717      do jcoeff=1,ncoeff_fit
1718        jcoeff_tmp = list_coeffs(jcoeff)
1719        enu = energy_coeffs(jcoeff_tmp,itime)
1720 !       etmpA =  emu*enu
1721 !       A(icoeff,jcoeff) = A(icoeff,jcoeff) + efact*etmpA
1722      end do
1723      etmpB = etmpB + energy_diff(itime)*emu / (sqomega(itime)**3)
1724      etmpB = zero ! REMOVE THIS LINE TO TAKE INTO ACOUNT THE ENERGY     
1725      
1726 !    Fill forces
1727      do ia=1,natom
1728        do mu=1,3
1729          fmu = fcart_coeffs_tmp(mu,ia,itime)
1730          do jcoeff=1,ncoeff_fit
1731            jcoeff_tmp = list_coeffs(jcoeff)
1732            fnu = fcart_coeffs(mu,ia,jcoeff_tmp,itime)
1733            ftmpA =  fmu*fnu
1734            A(icoeff,jcoeff) = A(icoeff,jcoeff) + ffact*ftmpA
1735          end do
1736          ftmpB = ftmpB + fcart_diff(mu,ia,itime)*fmu
1737        end do !End loop dir
1738      end do !End loop natom
1739 !    Fill stresses
1740      do mu=1,6
1741        smu = strten_coeffs(mu,itime,icoeff_tmp)
1742        do jcoeff=1,ncoeff_fit
1743          jcoeff_tmp = list_coeffs(jcoeff)
1744          snu = strten_coeffs(mu,itime,jcoeff_tmp)
1745          stmpA =  sqomega(itime)*smu*snu
1746          A(icoeff,jcoeff) = A(icoeff,jcoeff) + sfact*stmpA
1747        end do
1748        stmpB = stmpB + sqomega(itime)*strten_diff(mu,itime)*smu
1749      end do !End loop stress dir
1750    end do ! End loop time
1751    B(icoeff,1) = B(icoeff,1) + ffact*ftmpB + sfact*stmpB + efact*etmpB
1752  end do ! End loop icoeff
1753 
1754 !2-Solve Ax=B
1755 !OLD VERSION..
1756 ! call dgesvx(FACT,TRANS,N,NRHS,A,LDA,AF,LDAF,IPIV,EQUED,R,C,B,LDB,coefficients,LDX,&
1757 !             RCOND,FERR,BERR,WORK,IWORK,INFO)
1758 !U is nonsingular
1759 ! if (INFO==N+1) then
1760 !   coefficients = zero
1761 ! end if
1762  call DSGESV(N,NRHS,A,LDA,IPIV,B,LDB,coefficients,LDX,WORK,SWORK,ITER,INFO)
1763 
1764 !other routine
1765 ! call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
1766 ! coefficients = B(:,NRHS)
1767  !U is nonsingular
1768  if (INFO==N+2) then
1769    coefficients = zero
1770  end if
1771 
1772  if(any(abs(coefficients)>1.0E10))then
1773    INFO = 1
1774    coefficients = zero
1775  end if
1776 
1777  info_out = INFO
1778 
1779  ABI_DEALLOCATE(AF)
1780  ABI_DEALLOCATE(IPIV)
1781  ABI_DEALLOCATE(R)
1782  ABI_DEALLOCATE(C)
1783  ABI_DEALLOCATE(FERR)
1784  ABI_DEALLOCATE(BERR)
1785  ABI_DEALLOCATE(WORK)
1786  ABI_DEALLOCATE(IWORK)
1787  ABI_DEALLOCATE(SWORK)
1788  ABI_DEALLOCATE(A)
1789  ABI_DEALLOCATE(B)
1790 
1791 end subroutine fit_polynomial_coeff_solve

m_fit_polynomial_coeff/fit_polynomial_printSystemFiles [ Functions ]

[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]

NAME

 fit_polynomial_printSystemFiles

FUNCTION

 Print the files for the fitting script

INPUTS

 eff_pot<type(effective_potential)> = effective potential
 hist<type(abihist)> = datatype with the  history of the MD

OUTPUT

PARENTS

      multibinit

CHILDREN

      destroy_supercell,generelist,init_supercell,xred2xcart

SOURCE

2401 subroutine fit_polynomial_printSystemFiles(eff_pot,hist)
2402 
2403 
2404 !This section has been created automatically by the script Abilint (TD).
2405 !Do not modify the following lines by hand.
2406 #undef ABI_FUNC
2407 #define ABI_FUNC 'fit_polynomial_printSystemFiles'
2408 !End of the abilint section
2409 
2410  implicit none
2411 
2412 !Arguments ------------------------------------
2413 !scalars
2414 !arrays
2415  type(effective_potential_type), intent(in) :: eff_pot
2416  type(abihist),intent(in) :: hist
2417 !Local variables-------------------------------
2418 !scalar
2419  integer :: ia,ib,ib1,ii,jj,irpt,kk,ll,mu,nu,nstep,nshift
2420  integer :: natom_uc
2421  integer :: unit_born=22,unit_epsiloninf=23,unit_md=24
2422  integer :: unit_harmonic=25,unit_ref=26,unit_strain=27,unit_sym=28
2423 !arrays
2424  integer,allocatable :: typat_order(:),typat_order_uc(:)
2425  integer, dimension(3)  :: A,ncell
2426  real(dp), allocatable :: xcart(:,:),fcart(:,:)
2427  character(len=500) :: msg
2428  type(supercell_type) :: supercell
2429 ! *************************************************************************
2430 
2431 !Create new supercell corresponding to the MD
2432  ncell = (/2,2,2/)
2433  call init_supercell(eff_pot%crystal%natom, (/ncell(1),0,0,  0,ncell(2),0,  0,0,ncell(3)/),&
2434 &                    eff_pot%crystal%rprimd,eff_pot%crystal%typat,&
2435 &                    eff_pot%crystal%xcart,eff_pot%crystal%znucl, supercell)
2436 
2437 !allocation of array
2438  ABI_ALLOCATE(xcart,(3,supercell%natom))
2439  ABI_ALLOCATE(fcart,(3,supercell%natom))
2440  ABI_ALLOCATE(typat_order,(supercell%natom))
2441  ABI_ALLOCATE(typat_order_uc,(eff_pot%crystal%natom))
2442 
2443  A = (/ 2, 3, 1/)
2444 
2445  nshift = product(ncell)
2446  natom_uc = eff_pot%crystal%natom
2447 !Fill the typat_order array:
2448 !In the fit script the atom must be in the order 11111 222222 33333 ..
2449 !and the order of the atom can not be change in the fit script,
2450 !we transform into the format of the script
2451  ib = 1
2452  ib1= 1
2453  do ii=1,eff_pot%crystal%ntypat
2454    jj = A(ii)
2455    do kk=1,natom_uc
2456      if(supercell%typat(kk)==jj)then
2457        typat_order_uc(ib1) = kk
2458        ib1 = ib1 + 1
2459        do ll=1,nshift
2460          ia = (ll-1)*natom_uc + kk
2461          typat_order(ib) = ia
2462          ib = ib + 1
2463        end do
2464      end if
2465    end do
2466  end do
2467 
2468 ! BORN CHARGES FILE
2469  if (open_file('system/Born_Charges',msg,unit=unit_born,form="formatted",&
2470 &    status="replace",action="write") /= 0) then
2471    MSG_ERROR(msg)
2472  end if
2473  do ii=1,eff_pot%crystal%ntypat
2474    jj = A(ii)
2475    do ia=1,eff_pot%crystal%natom
2476      if(eff_pot%crystal%typat(ia)==jj)then
2477        write(unit_born,'(i2,a,1F10.5)') ia,"    ",eff_pot%crystal%amu(eff_pot%crystal%typat(ia))
2478        do mu=1,3
2479          WRITE(unit_born,'(a,3(F23.14))') "     ",eff_pot%harmonics_terms%zeff(:,mu,ia)
2480        end do
2481      end if
2482    end do
2483  end do
2484 
2485 !DIELECTRIC TENSOR FILE
2486  if (open_file('system/Dielectric_Tensor',msg,unit=unit_epsiloninf,form="formatted",&
2487 &    status="replace",action="write") /= 0) then
2488    MSG_ERROR(msg)
2489  end if
2490  do mu=1,3
2491    WRITE(unit_epsiloninf,'(3(F23.14))') eff_pot%harmonics_terms%epsilon_inf(:,mu)
2492  end do
2493 
2494 
2495 !REFERENCE STRUCTURE FILE
2496  if (open_file('system/Reference_structure',msg,unit=unit_ref,form="formatted",&
2497 &    status="replace",action="write") /= 0) then
2498    MSG_ERROR(msg)
2499  end if
2500 
2501  write(unit_ref,'("Energy (Hartree)")')
2502  write(unit_ref,'("================")')
2503  write(unit_ref,'(F23.14)') (hist%etot(1)/nshift)
2504  write(unit_ref,'("")')
2505  write(unit_ref,'("Cell vectors")')
2506  write(unit_ref,'("============")')
2507  do jj=1,3
2508    write(unit_ref,'(3(F22.14))') (supercell%rprimd(:,jj))
2509  end do
2510 
2511  write(unit_ref,'("")')
2512  write(unit_ref,'("Atomic positions (Bohr radius)")')
2513  write(unit_ref,'("==============================")')
2514 
2515  do ia=1,supercell%natom
2516    write(unit_ref,'(3(F23.14))') supercell%xcart(:,typat_order(ia))
2517  end do
2518 
2519 !Harmonic XML file
2520  if (open_file('system/harmonic.xml',msg,unit=unit_harmonic,form="formatted",&
2521 &     status="replace",action="write") /= 0) then
2522    MSG_ERROR(msg)
2523  end if
2524 
2525 !Write header
2526  write(unit_harmonic,'("<?xml version=""1.0"" ?>")')
2527  write(unit_harmonic,'("<name>")')
2528 
2529  do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt
2530    if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol9)) then
2531      write(unit_harmonic,'("  <local_force_constant units=""hartree/bohrradius**2"">")')
2532      write(unit_harmonic,'("    <data>")')
2533      do ia=1,eff_pot%crystal%natom
2534        do mu=1,3
2535          do ib=1,eff_pot%crystal%natom
2536            do  nu=1,3
2537              write(unit_harmonic,'(F22.14)', advance="no")&
2538 &                 (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,typat_order_uc(ia),&
2539 &                                                              nu,typat_order_uc(ib),irpt))
2540            end do
2541          end do
2542          write(unit_harmonic,'(a)')''
2543        end do
2544      end do
2545      write(unit_harmonic,'("    </data>")')
2546      write(unit_harmonic,'("    <cell>")')
2547      write(unit_harmonic,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
2548      write(unit_harmonic,'("    </cell>")')
2549      write(unit_harmonic,'("  </local_force_constant>")')
2550    end if
2551  end do
2552  write(unit_harmonic,'("</name>")')
2553 
2554 !STRAIN FILE
2555  if (open_file('system/Strain_Tensor',msg,unit=unit_strain,form="formatted",&
2556 &     status="replace",action="write") /= 0) then
2557    MSG_ERROR(msg)
2558  end if
2559  write(unit_strain,'(6(F23.14))') (eff_pot%harmonics_terms%elastic_constants)
2560 
2561 ! SYM FILE
2562  if (open_file('system/symmetry_operations',msg,unit=unit_sym,form="formatted",&
2563 &     status="replace",action="write") /= 0) then
2564    MSG_ERROR(msg)
2565  end if
2566  write(unit_sym,'("(x,y,z)  (y,-x,z) (z,x,y) (y,z,x) (x,z,y) (y,x,z) (z,y,x) (x,-y,-z) (z,-x,-y)",&
2567 &                " (y,-z,-x) (x,-z,-y) (y,-x,-z) (z,-y,-x) (-x,y,-z) (-z,x,-y) (-y,z,-x) (-x,z,-y)",&
2568 &                " (-y,x,-z) (-z,y,-x) (-x,-y,z) (-z,-x,y) (-y,-z,x) (-x,-z,y) (-y,-x,z) (-z,-y,x)",&
2569 &                " (-x,-y,-z) (-z,-x,-y) (-y,-z,-x) (-x,-z,-y) (-y,-x,-z) (-z,-y,-x) (-x,y,z)",&
2570 &                " (-z,x,y) (-y,z,x) (-x,z,y) (-y,x,z) (-z,y,x) (x,-y,z) (z,-x,y) (y,-z,x) (x,-z,y)",&
2571 &                " (z,-y,x) (x,y,-z) (z,x,-y) (y,z,-x) (x,z,-y) (y,x,-z) (z,y,-x)")')
2572 
2573 
2574 !MD file
2575  nstep = hist%mxhist
2576  if (open_file('system/Molecular_dynamic',msg,unit=unit_md,form="formatted",&
2577 &     status="replace",action="write") /= 0) then
2578    MSG_ERROR(msg)
2579  end if
2580  do ii=1,nstep
2581    write(unit_md,'(I5)') ii-1
2582    write(unit_md,'(F22.14)') hist%etot(ii)/nshift
2583    do jj=1,3
2584      write(unit_md,'(3(F22.14))') (hist%rprimd(:,jj,ii))
2585    end do
2586 !  Set xcart and fcart for this step
2587    call xred2xcart(supercell%natom,hist%rprimd(:,:,ii),&
2588 &                  xcart,hist%xred(:,:,ii))
2589 
2590    fcart(:,:) = hist%fcart(:,:,ii)
2591 
2592    do ia=1,supercell%natom
2593      write(unit_md,'(3(E22.14),3(E22.14))') xcart(:,typat_order(ia)),fcart(:,typat_order(ia))
2594    end do
2595    write(unit_md,'(6(E22.14))') hist%strten(:,ii)
2596  end do
2597 
2598 !Close files
2599  close(unit_ref)
2600  close(unit_born)
2601  close(unit_harmonic)
2602  close(unit_epsiloninf)
2603  close(unit_md)
2604  close(unit_strain)
2605  close(unit_sym)
2606 
2607 !Deallocation array
2608  ABI_DEALLOCATE(typat_order)
2609  ABI_DEALLOCATE(typat_order_uc)
2610  ABI_DEALLOCATE(xcart)
2611  ABI_DEALLOCATE(fcart)
2612  call destroy_supercell(supercell)
2613 
2614 end subroutine fit_polynomial_printSystemFiles