TABLE OF CONTENTS


ABINIT/m_mover_effpot [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mover_effpot

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_mover_effpot
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  implicit none
34 
35  private

ABINIT/mover_effpot [ Functions ]

[ Top ] [ Functions ]

NAME

 mover_effpot

FUNCTION

 this routine is driver for using mover with effective potential

COPYRIGHT

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

INPUTS

  inp = input of multibinit
  effective_potential =  effective potential of the reference structure
  option = flag for the option:
                         -1  = Bound the anharmonic part
                          12 = NVT simulation
                          13 = NPT simulation

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      multibinit

CHILDREN

      alloc_copy,destroy_mpi_enreg,destroy_results_gs,dtset_free
      effective_potential_setcoeffs,effective_potential_setsupercell
      fit_polynomial_coeff_fit,fit_polynomial_coeff_getpositive,generelist
      init_results_gs,mover,polynomial_coeff_free,polynomial_coeff_getnorder
      polynomial_coeff_init,polynomial_coeff_setcoefficient
      polynomial_coeff_writexml,scfcv_destroy,wrtout,xcart2xred,xmpi_barrier
      xred2xcart

SOURCE

  86 subroutine mover_effpot(inp,filnam,effective_potential,option,comm,hist)
  87 
  88  use defs_basis
  89  use defs_abitypes
  90  use m_abicore
  91  use defs_datatypes
  92  use m_errors
  93  use m_abimover
  94  use m_build_info
  95  use m_scf_history
  96  use defs_wvltypes
  97  use m_xmpi
  98  use m_abimover
  99  use m_phonons
 100  use m_strain
 101  use m_effective_potential_file
 102  use m_supercell
 103  use m_psps
 104  use m_args_gs
 105 
 106  use m_geometry, only : xcart2xred, xred2xcart
 107  use m_multibinit_dataset, only : multibinit_dtset_type
 108  use m_effective_potential,only : effective_potential_type
 109  use m_fit_polynomial_coeff, only : polynomial_coeff_writeXML
 110  use m_fit_polynomial_coeff, only : fit_polynomial_coeff_fit,genereList
 111  use m_fit_polynomial_coeff, only : fit_polynomial_coeff_getPositive,fit_polynomial_coeff_getCoeffBound
 112  use m_electronpositron,   only : electronpositron_type
 113  use m_polynomial_coeff,only : polynomial_coeff_getNorder
 114 ! use m_pawang,       only : pawang_type, pawang_free
 115 ! use m_pawrad,       only : pawrad_type, pawrad_free
 116 ! use m_pawtab,       only : pawtab_type, pawtab_nullify, pawtab_free
 117 ! use m_pawxmlps, only : paw_setup, ipsp2xml, rdpawpsxml, &
 118 !&                       paw_setup_copy, paw_setup_free, getecutfromxml
 119  use m_dtset,  only : dtset_free
 120  use m_abihist, only : abihist
 121  use m_ifc
 122  use m_ewald
 123  use m_mpinfo,           only : init_mpi_enreg,destroy_mpi_enreg
 124  use m_copy            , only : alloc_copy
 125  use m_electronpositron, only : electronpositron_type
 126  use m_scfcv,            only : scfcv_t, scfcv_run,scfcv_destroy
 127  use m_results_gs,       only : results_gs_type,init_results_gs,destroy_results_gs
 128  use m_mover,            only : mover
 129 
 130 !This section has been created automatically by the script Abilint (TD).
 131 !Do not modify the following lines by hand.
 132 #undef ABI_FUNC
 133 #define ABI_FUNC 'mover_effpot'
 134 !End of the abilint section
 135 
 136 implicit none
 137 
 138 !Arguments --------------------------------
 139 !scalar
 140  integer, intent(in) :: option,comm
 141 !array
 142  type(multibinit_dtset_type),intent(in) :: inp
 143  type(effective_potential_type),intent(inout)  :: effective_potential
 144  character(len=fnlen),intent(in) :: filnam(15)
 145  type(abihist),optional,intent(inout):: hist
 146 !Local variables-------------------------------
 147 !scalar
 148  integer :: filetype,icoeff_bound,ii
 149 !integer :: iexit,initialized
 150  integer :: jj,kk,nproc,ncoeff,nmodels,ncoeff_bound,ncoeff_bound_tot,ncoeff_max
 151  integer :: model_bound,model_ncoeffbound,my_rank
 152 !integer :: mtypalch,,npsp,paw_size,type
 153 !integer,save :: paw_size_old=-1
 154  real(dp):: cutoff,freq_q,freq_b,qmass,bmass
 155 ! real(dp):: time_q,time_b
 156  logical :: iam_master
 157  integer, parameter:: master=0
 158  logical :: verbose,writeHIST
 159 !real(dp) :: cpui
 160 !character(len=6) :: codvsn
 161 
 162 !TEST_AM
 163 ! integer :: ia,mu,rand_seed = 5
 164 ! real(dp):: mass_ia,rescale_vel,sum_mass,v2gauss
 165 !TEST_AM
 166 ! Set array dimensions
 167  character(len=500) :: message
 168  type(MPI_type),target :: mpi_enreg
 169  type(dataset_type),target :: dtset
 170  type(scfcv_t) :: scfcv_args
 171  type(datafiles_type),target :: dtfil
 172  integer,target :: zero_integer
 173  type(ab_xfh_type) :: ab_xfh
 174  type(results_gs_type),target :: results_gs
 175  type(pseudopotential_type),target :: psps
 176 !arrays
 177 !no_abirules
 178  integer :: sc_size(3),sc_size_TS(3)
 179  integer,pointer :: indsym(:,:,:)
 180  integer,allocatable :: listcoeff(:),listcoeff_bound(:,:),list_tmp(:),list_bound(:,:)
 181  integer,allocatable :: isPositive(:)
 182  integer,allocatable :: symrel(:,:,:)
 183 !integer,allocatable :: npwtot(:)
 184  real(dp) :: acell(3)
 185 !real(dp) :: ecut_tmp(3,2,10)
 186  real(dp),allocatable :: coeff_values(:,:)
 187  real(dp),pointer :: rhog(:,:),rhor(:,:)
 188  real(dp),allocatable :: tnons(:,:)
 189  real(dp),allocatable :: xred(:,:),xred_old(:,:),xcart(:,:)
 190  real(dp),allocatable :: fred(:,:),fcart(:,:)
 191  real(dp),allocatable :: vel(:,:)
 192  real(dp) :: vel_cell(3,3),rprimd(3,3)
 193  type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_all,coeffs_tmp,coeffs_bound
 194  character(len=fnlen) :: filename
 195 !character(len=fnlen) :: filename_psp(3)
 196  type(electronpositron_type),pointer :: electronpositron
 197 ! type(pspheader_type),allocatable :: pspheads(:)
 198 ! type(pawrad_type),allocatable :: pawrad(:)
 199 ! type(pawtab_type),allocatable :: pawtab(:)
 200 ! type(args_gs_type) :: args_gs
 201 !type(wvl_data) :: wvl
 202 !type(pawang_type) :: pawang
 203 !type(scf_history_type) :: scf_history
 204 
 205 !******************************************************************
 206 
 207 !MPI variables
 208  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
 209  iam_master = (my_rank == master)
 210 
 211  write(message, '(a,(80a),a)') ch10,&
 212 & ('=',ii=1,80),ch10
 213  call wrtout(ab_out,message,'COLL')
 214  call wrtout(std_out,message,'COLL')
 215 
 216 !*******************************************************************
 217 ! 1 Generate and check supercell for the dynamics
 218 !*******************************************************************
 219 
 220 !a new supercell is compute
 221 !Initialisaton of variable
 222  if(option == -1.or.option == -2) then
 223 !  Bound process option
 224    sc_size(:) = inp%bound_cell
 225  else if(option == -3) then
 226 !  Heff option
 227    sc_size(:) = (/1,1,1/)
 228  else
 229 !  Normal dynamics
 230    sc_size(:) = inp%ncell
 231  end if
 232 
 233  if(option/=0)then
 234 
 235 !  acell is always set to one, only rprimd is used for the effective potential
 236    acell = one
 237    rprimd = effective_potential%crystal%rprimd
 238 
 239    ABI_ALLOCATE(xred,(3,effective_potential%crystal%natom))
 240    ABI_ALLOCATE(xcart,(3,effective_potential%crystal%natom))
 241 
 242 !  convert new xcart
 243    call xcart2xred(effective_potential%crystal%natom,effective_potential%crystal%rprimd,&
 244 &   effective_potential%crystal%xcart,xred)
 245    call xred2xcart(effective_potential%crystal%natom, rprimd, xcart, xred)
 246 !  Generate supercell for the simulation
 247    call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
 248 
 249    ABI_DEALLOCATE(xred)
 250    ABI_DEALLOCATE(xcart)
 251 
 252 !***************************************************************
 253 !1 Convert some parameters into the structures used by mover.F90
 254 !***************************************************************
 255 !NOTE:ARGUMENTS OF MOVER SHOULD BE CLEAN
 256 !     We may just need to provide AB_MOVER wich is the main object
 257 !     for mover and set  scfcv_args as an optional and depending on
 258 !     the kind of calculation (abinit or multibinit), we provide
 259 !     to mover scfcv_args or effective_potential...
 260 !***************************************************************
 261 !  Free dtset
 262    call dtset_free(dtset)
 263 
 264 !  Set mpi_eng
 265    mpi_enreg%comm_cell  = comm
 266    mpi_enreg%me = my_rank
 267 
 268 !  Set the abinit dataset for mover with fake values
 269 !  Scalar
 270    dtset%dmft_entropy = 0
 271    dtset%nctime = inp%nctime ! NetCdf TIME between output of molecular dynamics informations
 272    dtset%delayperm = 0  ! DELAY between trials to PERMUTE atoms
 273    dtset%dilatmx = one  ! DILATation : MaXimal value
 274    dtset%chkdilatmx = 0 ! No check on dilatmx is needed in multibilint
 275    dtset%diismemory = 8 ! Direct Inversion in the Iterative Subspace MEMORY
 276    dtset%friction = 0.0001d0 ! internal FRICTION coefficient
 277    dtset%goprecon = 0   ! Geometry Optimization PREconditioner equations
 278    dtset%istatr = 0     ! Integer for STATus file SHiFT
 279    dtset%jellslab = 0   ! include a JELLium SLAB in the cell
 280    dtset%mqgrid = 0     ! Maximum number of Q-space GRID points for pseudopotentials
 281    dtset%mqgriddg = 0   ! Maximum number of Q-wavevectors for the 1-dimensional GRID
 282                         ! for the Double Grid in PAW
 283    dtset%mdwall = 10000d0 ! Molecular Dynamics WALL location
 284    dtset%ntypalch = 0   ! Number of TYPe of atoms that are "ALCHemical"
 285    dtset%natom = effective_potential%supercell%natom
 286    dtset%ntypat = effective_potential%crystal%ntypat
 287    dtset%npspalch = effective_potential%crystal%ntypat
 288    dtset%nconeq = 0     ! Number of CONstraint EQuations
 289    dtset%noseinert = 1.d-5 ! NOSE INERTia factor
 290    dtset%nnos = inp%nnos   ! Number of nose masses Characteristic
 291    do ii=1,3 ! Only copy the diagonal part
 292              ! might be adapt with ngqpt(9) instead of ngqpt(3) because is needed in wght9.f
 293      dtset%ph_ngqpt(ii) = inp%ngqpt((ii))
 294    end do
 295    dtset%ph_nqshift = inp%nqshft
 296    dtset%nsym = 1       ! Number of SYMmetry operations
 297    dtset%prtxml = 0     ! print the xml
 298    dtset%signperm = 1   ! SIGN of PERMutation potential
 299    dtset%strprecon = 1  ! STRess PRECONditioner
 300    dtset%supercell_latt(:,:) = 0
 301    do ii=1,3
 302      dtset%supercell_latt(ii,ii) = sc_size(ii)
 303    end do
 304    dtset%tolmxf = 2.0d-5
 305    dtset%tsmear = 0.009500446 !
 306    dtset%vis = 100      ! VIScosity
 307    dtset%usewvl = 0     !
 308    dtset%useylm = 0     !
 309 
 310 !    if(option == -3) then
 311 !      write(message,'(a)')' Read the DDB file to fill the dtset array'
 312 !      call wrtout(std_out,message,"COLL")
 313 ! !    Copy real informtions from the ddb
 314 !      call effective_potential_file_getType(filnam(3),type)
 315 !      if(type /= 1) then
 316 !        write(message, '(5a)' )&
 317 ! &         ' You need to provide DDB file in the input to compute ahnarmonic',ch10,&
 318 ! &         ' part of effective Hamiltionian',ch10,&
 319 ! &         'Action: add DDB file in the inputs'
 320 !        MSG_BUG(message)
 321 !      end if
 322 !      call ddb_to_dtset(comm, dtset,filnam(3),psps)
 323 !      ABI_ALLOCATE(dtset%kptns,(3,dtset%nkpt))
 324 !      dtset%kptns(:,:) = dtset%kpt(:,:)
 325 !      ABI_ALLOCATE(dtset%istwfk,(dtset%nkpt))
 326 !      dtset%istwfk(:) = 1
 327 
 328 !    else
 329 !    Need to init some values
 330    ABI_ALLOCATE(symrel,(3,3,dtset%nsym))
 331    symrel = 1
 332    call alloc_copy(symrel,dtset%symrel)
 333    ABI_ALLOCATE(tnons,(3,dtset%nsym))
 334    tnons = zero
 335    call alloc_copy(tnons,dtset%tnons)
 336    call alloc_copy(effective_potential%supercell%typat,dtset%typat)
 337    call alloc_copy(effective_potential%crystal%znucl,dtset%znucl)
 338    ABI_DEALLOCATE(symrel)
 339    ABI_DEALLOCATE(tnons)
 340 !   end if
 341 
 342    !array
 343    ABI_ALLOCATE(dtset%iatfix,(3,dtset%natom)) ! Indices of AToms that are FIXed
 344    dtset%iatfix = 0
 345    dtset%goprecprm(:) = zero !Geometry Optimization PREconditioner PaRaMeters equations
 346    ABI_ALLOCATE(dtset%prtatlist,(dtset%natom)) !PRinT by ATom LIST of ATom
 347    dtset%prtatlist(:) = 0
 348    ABI_ALLOCATE(dtset%mixalch_orig,(dtset%npspalch,dtset%ntypalch,1))
 349    dtset%mixalch_orig(:,:,:)=zero
 350    ABI_ALLOCATE(dtset%ph_qshift,(3,dtset%ph_nqshift))
 351    dtset%ph_qshift = inp%q1shft
 352    if(option  > 0)then
 353      verbose = .TRUE.
 354      writeHIST = .TRUE.
 355      dtset%dtion = inp%dtion  ! Delta Time for IONs
 356      dtset%ionmov = inp%dynamics  ! Number for the dynamic
 357      dtset%ntime = inp%ntime  ! Number of TIME steps
 358      dtset%optcell = inp%optcell    ! OPTimize the CELL shape and dimensions Characteristic
 359      dtset%restartxf = inp%restartxf  ! RESTART from (X,F) history
 360      dtset%mdtemp(1) = inp%temperature   !Molecular Dynamics Temperatures
 361      dtset%mdtemp(2) = inp%temperature   !Molecular Dynamics Temperatures
 362      dtset%strtarget(1:6) = -1 * inp%strtarget(1:6) / 29421.033d0 ! STRess TARGET
 363    else if(option == -1.or.option == -2) then
 364 !    Set default for the fit
 365      verbose = .false.
 366      writeHIST = .false.
 367      dtset%restartxf = 0  ! RESTART from (X,F) history
 368      dtset%dtion = inp%dtion  ! Delta Time for IONs
 369      dtset%ionmov = 13  ! Number for the dynamic
 370      dtset%ntime = inp%bound_step  ! Number of TIME steps
 371      dtset%optcell = 2    ! OPTimize the CELL shape and dimensions Characteristic
 372      dtset%mdtemp(1) = inp%bound_temp   !Molecular Dynamics Temperatures
 373      dtset%mdtemp(2) = inp%bound_temp !Molecular Dynamics Temperatures
 374      dtset%strtarget(1:6) = zero
 375    end if
 376 
 377 !  Set the barostat and thermonstat if ionmov == 13
 378    if(dtset%ionmov == 13)then
 379 
 380 !    Select frequency of the barostat as a function of temperature
 381 !    For small temperature, we need huge barostat and inversely
 382      ! if(dtset%mdtemp(1) <= 10) then
 383      !   freq_q = 0.002
 384      !   freq_b = 0.0002
 385      ! else  if(dtset%mdtemp(1) <= 50) then
 386      !   freq_q = 0.02
 387      !   freq_b = 0.002
 388      ! else  if(dtset%mdtemp(1) > 50.and.dtset%mdtemp(1) < 300) then
 389      !   freq_q = 0.1
 390      !   freq_b = 0.01
 391      ! else
 392      !   freq_q = 0.2
 393      !   freq_b = 0.02
 394      ! end if
 395 
 396 !TEST_AM_old way
 397 !     freq_q = 0.1
 398 !     freq_b = 0.01
 399 !     qmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_q**2)
 400 !     bmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_b**2)
 401 !TEST_AM
 402 
 403 
 404 !TEST_AM
 405      freq_q = 800  / Ha_cmm1
 406      freq_b = 800  / Ha_cmm1
 407      qmass = 10 * dtset%natom * kb_HaK * dtset%mdtemp(1) / (freq_q**2)
 408      bmass = 10000*qmass
 409 !TEST_AM
 410 
 411 
 412      if(dtset%nnos==0) then
 413        dtset%nnos = 1
 414        ABI_ALLOCATE(dtset%qmass,(dtset%nnos))
 415        dtset%qmass(:)  = qmass
 416        write(message,'(3a,F30.10,a)')&
 417 &       ' WARNING: nnos is set to zero in the input',ch10,&
 418 &       '          value by default for qmass: ',dtset%qmass(:),ch10
 419        if(verbose)call wrtout(std_out,message,"COLL")
 420      else
 421        ABI_ALLOCATE(dtset%qmass,(dtset%nnos)) ! Q thermostat mass
 422        dtset%qmass(:) = inp%qmass(:)
 423      end if
 424      if (abs(inp%bmass) < tol10) then
 425        dtset%bmass = bmass
 426        write(message,'(3a,F30.10,a)')&
 427 &       ' WARNING: bmass is set to zero in the input',ch10,&
 428 &       '          value by default for bmass: ',dtset%bmass,ch10
 429        if(verbose)call wrtout(std_out,message,"COLL")
 430      else
 431        dtset%bmass = inp%bmass  ! Barostat mass
 432      end if
 433    end if
 434 
 435 !  Do a check
 436    if(dtset%ionmov == 27)then
 437      call effective_potential_file_getType(filnam(3),filetype)
 438      if(filetype /= 1)then
 439        write(message, '(5a)' )&
 440 &       ' The file ',trim(filnam(3)),' is not a DDB',ch10,&
 441 &       ' It is not compatible with ionmov 27'
 442        MSG_ERROR(message)
 443      end if
 444 
 445    end if
 446 
 447 !  set psps
 448    psps%useylm = dtset%useylm
 449 !    if(option == -3)then
 450 !      mtypalch = 0
 451 !      npsp = dtset%ntypat
 452 !      call psps_free(psps)
 453 !      filename_psp(1) = "/home/alex/calcul/psp/Sr.LDA_PW-JTH.xml"
 454 !      filename_psp(2) = "/home/alex/calcul/psp/Ti.LDA_PW-JTH.xml"
 455 !      filename_psp(3) = "/home/alex/calcul/psp/O.LDA_PW-JTH.xml"
 456 !      ABI_DATATYPE_ALLOCATE(pspheads,(npsp))
 457 !      call inpspheads(filename_psp,npsp,pspheads,ecut_tmp)
 458 !      call psps_init_global(mtypalch, npsp, psps, pspheads)
 459 !      call psps_init_from_dtset(dtset, 1, psps, pspheads)
 460 !    end if
 461 
 462 ! !  The correct dimension of pawrad/tab is ntypat. In case of alchemical psps
 463 ! !  pawrad/tab(ipsp) is invoked with ipsp<=npsp. So, in order to avoid any problem,
 464 ! !  declare pawrad/tab at paw_size=max(ntypat,npsp).
 465 !    paw_size=0;if (psps%usepaw==1) paw_size=max(dtset%ntypat,dtset%npsp)
 466 !    if (paw_size/=paw_size_old) then
 467 !      if (paw_size_old/=-1) then
 468 !        call pawrad_free(pawrad)
 469 !        call pawtab_free(pawtab)
 470 !        ABI_DATATYPE_DEALLOCATE(pawrad)
 471 !        ABI_DATATYPE_DEALLOCATE(pawtab)
 472 !      end if
 473 !      ABI_DATATYPE_ALLOCATE(pawrad,(paw_size))
 474 !      ABI_DATATYPE_ALLOCATE(pawtab,(paw_size))
 475 !      call pawtab_nullify(pawtab)
 476 !      paw_size_old=paw_size
 477 !    end if
 478 
 479 !  set args_gs
 480 !    if (option == -3 then
 481 !      call args_gs_init(args_gs, &
 482 ! &       effective_potential%crystal%amu(:),dtset%mixalch_orig(:,:,1),&
 483 ! &       dtset%dmatpawu(:,:,:,:,1),dtset%upawu(:,1),dtset%jpawu(:,1),&
 484 ! &       dtset%rprimd_orig(:,:,1))
 485 !      ABI_ALLOCATE(npwtot,(dtset%nkpt))
 486 !    end if
 487 !  initialisation of results_gs
 488    call init_results_gs(dtset%natom,1,results_gs)
 489 
 490 !  Set the pointers of scfcv_args
 491    zero_integer = 0
 492    scfcv_args%dtset     => dtset
 493    ABI_ALLOCATE(indsym,(4,dtset%nsym,dtset%natom))
 494    indsym = 0
 495    scfcv_args%indsym => indsym
 496    scfcv_args%mpi_enreg => mpi_enreg
 497    scfcv_args%ndtpawuj  => zero_integer
 498    scfcv_args%results_gs => results_gs
 499    scfcv_args%psps => psps
 500 !  Set other arguments of the mover.F90 routines
 501 
 502 !  Set the dffil structure
 503    dtfil%filnam_ds(1:2)=filnam(1:2)
 504    dtfil%filnam_ds(3)=""
 505    dtfil%filnam_ds(4)=filnam(2)
 506    dtfil%filstat='_STATUS'
 507    nullify (electronpositron)
 508    ABI_ALLOCATE(rhog,(2,1))
 509    ABI_ALLOCATE(rhor,(2,1))
 510 
 511 !  Initialize xf history (should be put in inwffil)
 512 !  Not yet implemented for ionmov 2 3 10 11 22 (memory problem...)
 513 !  ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5
 514    ab_xfh%nxfh = 0
 515    ab_xfh%mxfh = 1
 516    ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh))
 517    if (any((/2,3,10,11,22/)==dtset%ionmov)) then
 518      write(message, '(3a)' )&
 519 &     ' This dynamics can not be used with effective potential',ch10,&
 520 &     'Action: correct dynamics input'
 521      MSG_BUG(message)
 522    end if
 523 
 524 !***************************************************************
 525 !2  initialization of the structure for the dynamics
 526 !***************************************************************
 527 
 528    if (allocated(dtset%rprimd_orig)) then
 529      ABI_DEALLOCATE(dtset%rprimd_orig)
 530    end if
 531    ABI_ALLOCATE(dtset%rprimd_orig,(3,3,1))
 532    dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
 533 
 534 
 535    ABI_ALLOCATE(xred,(3,dtset%natom))
 536    ABI_ALLOCATE(xred_old,(3,dtset%natom))
 537    ABI_ALLOCATE(vel,(3,dtset%natom))
 538    ABI_ALLOCATE(fred,(3,dtset%natom))
 539    ABI_ALLOCATE(fcart,(3,dtset%natom))
 540 
 541    call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
 542 &   effective_potential%supercell%xcart,xred)
 543 
 544    xred_old = xred
 545    vel_cell(:,:) = zero
 546    vel(:,:)      = zero
 547 
 548 !*********************************************************
 549 !4   Call main routine for the bound process,
 550 !    monte carlo / molecular dynamics / project
 551 !*********************************************************
 552    if(option > 0)then
 553      !*************************************************************
 554      !  call mover in case of NPT or NVT simulation
 555      !*************************************************************
 556      write(message, '((80a),3a)' ) ('-',ii=1,80), ch10,&
 557 &     '-Monte Carlo / Molecular Dynamics ',ch10
 558      call wrtout(ab_out,message,'COLL')
 559      call wrtout(std_out,message,'COLL')
 560      call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
 561 &     rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
 562 &     effective_potential=effective_potential,filename_ddb=filnam(3),&
 563 &     verbose=verbose,writeHIST=writeHIST)
 564 
 565    else if(option== -1.or.option==-2)then
 566      !*************************************************************
 567      !   Try to bound the model
 568      !*************************************************************
 569      write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,&
 570 &     ' Try to bound the model',ch10,' Check if the model is bounded or not'
 571      call wrtout(ab_out,message,'COLL')
 572      call wrtout(std_out,message,'COLL')
 573 
 574 !    Try the model
 575      call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
 576 &     rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
 577 &     effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST)
 578 
 579      write(message, '(a)' ) ' => The model'
 580      if(effective_potential%anharmonics_terms%bounded)then
 581        write(message, '(2a)' ) trim(message),' is bound'
 582        call wrtout(std_out,message,'COLL')
 583        call wrtout(ab_out,message,'COLL')
 584      else
 585        write(message, '(2a)' ) trim(message),' is not bound'
 586        call wrtout(std_out,message,'COLL')
 587 
 588 
 589        if(option==-2)then
 590 
 591 !        Fill the list for the fixcoeff input of the fit_polynomial_coeff_fit routine
 592 !        Store the number of coefficients before adding other coeff for the bounding
 593          ncoeff = effective_potential%anharmonics_terms%ncoeff
 594          ABI_ALLOCATE(listcoeff,(ncoeff))
 595          do ii=1,ncoeff
 596            listcoeff(ii)=ii
 597          end do
 598 
 599          write(message, '(2a)')ch10,' Generate the list of addionnal terms with the fit process...'
 600          call wrtout(std_out,message,'COLL')
 601 
 602 !        Get the additional coeff
 603          call fit_polynomial_coeff_fit(effective_potential,(/0/),listcoeff,hist,1,&
 604 &         inp%bound_rangePower,0,inp%bound_maxCoeff,ncoeff,1,comm,cutoff_in=inp%bound_cutoff,&
 605 &         max_power_strain=2,verbose=.true.,positive=.true.,spcoupling=inp%bound_SPCoupling==1,&
 606 &         anharmstr=inp%bound_anhaStrain==1,only_even_power=.true.)
 607 
 608 !        Store the max number of coefficients after the fit process
 609          ncoeff_max = effective_potential%anharmonics_terms%ncoeff
 610 !        Store all the coefficients in coeffs_all
 611          ABI_DATATYPE_ALLOCATE(coeffs_all,(ncoeff_max))
 612          ABI_DATATYPE_ALLOCATE(coeffs_tmp,(ncoeff_max))
 613          do ii=1,ncoeff_max
 614            call polynomial_coeff_init(&
 615 &           effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
 616 &           effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
 617 &           coeffs_all(ii),&
 618 &           effective_potential%anharmonics_terms%coefficients(ii)%terms,&
 619 &           effective_potential%anharmonics_terms%coefficients(ii)%name,&
 620 &           check=.false.)
 621          end do
 622 
 623          model_ncoeffbound = 0
 624          do ii=1,ncoeff_max-ncoeff
 625 
 626            write(message, '(5a,I0,a)')ch10,'--',ch10,' Try to bound the model ',&
 627 &           'with ', ii,' additional term'
 628            if(ii>1)write(message,'(2a)') trim(message),'s'
 629            call wrtout(std_out,message,'COLL')
 630            call wrtout(ab_out,message,'COLL')
 631 
 632 !          Copy the new model in coeffs_tmp(jj)
 633 !          Free the coeffs_tmp array before
 634            do jj=1,ncoeff_max
 635              call polynomial_coeff_free(coeffs_tmp(jj))
 636            end do
 637 
 638            do jj=1,ncoeff+ii
 639              call polynomial_coeff_init(&
 640 &             coeffs_all(jj)%coefficient,&
 641 &             coeffs_all(jj)%nterm,&
 642 &             coeffs_tmp(jj),&
 643 &             coeffs_all(jj)%terms,&
 644 &             coeffs_all(jj)%name,&
 645 &             check=.false.)
 646            end do
 647 
 648            model_ncoeffbound = ii
 649 
 650 !          Reset the simulation and set the coefficients of the model
 651            call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,ncoeff+ii)
 652            call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,&
 653 &           -1,1,comm,verbose=.true.,positive=.false.)
 654            call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
 655            dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
 656            acell(:) = one
 657            call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
 658 &           effective_potential%supercell%xcart,xred)
 659            xred_old = xred
 660            vel_cell(:,:) = zero
 661            vel(:,:)      = zero
 662            fred(:,:)     = zero
 663            fcart(:,:)    = zero
 664 
 665 !          Run mover to check if the model is bound
 666            call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
 667 &           rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
 668 &           effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST)
 669            if(.not.effective_potential%anharmonics_terms%bounded)then
 670              write(message, '(2a)' ) ' => The model is not bounded'
 671            else
 672              write(message, '(2a)' ) ' => The model is bounded'
 673            end if
 674            call wrtout(std_out,message,'COLL')
 675 
 676 !          Exit if the model is bounded
 677            if(effective_potential%anharmonics_terms%bounded)  exit
 678 
 679          end do
 680 
 681        else
 682 !      Get the list of possible coefficients to bound the model
 683          cutoff = zero
 684          do ii=1,3
 685            cutoff = cutoff + effective_potential%crystal%rprimd(ii,ii)
 686          end do
 687          cutoff = cutoff / 3.0
 688 
 689 !        call fit_polynomial_coeff_getCoeffBound(effective_potential,coeffs_bound,&
 690 !&                                              hist,ncoeff_bound,comm,verbose=.true.)
 691 
 692 !TEST_AM!
 693          sc_size_TS = (/2,2,2/)
 694          call polynomial_coeff_getNorder(coeffs_bound,effective_potential%crystal,cutoff,&
 695 &         ncoeff_bound,ncoeff_bound_tot,inp%bound_rangePower,inp%bound_rangePower(2),2,sc_size_TS,&
 696 &         comm,anharmstr=inp%bound_anhaStrain==1,&
 697 &         spcoupling=inp%bound_SPCoupling==1,verbose=.false.,distributed=.false.,&
 698 &         only_even_power=.true.,only_odd_power=.false.)
 699 
 700          if(iam_master)then
 701            filename=trim(filnam(2))//"_boundcoeff.xml"
 702            call polynomial_coeff_writeXML(coeffs_bound,ncoeff_bound,filename=filename,newfile=.true.)
 703          end if
 704 !      wait
 705          call xmpi_barrier(comm)
 706 !      Store all the initial coefficients
 707          ncoeff = effective_potential%anharmonics_terms%ncoeff
 708          ABI_DATATYPE_ALLOCATE(coeffs_all,(ncoeff+ncoeff_bound))
 709          do ii=1,ncoeff
 710            call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
 711 &           effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
 712 &           coeffs_all(ii),&
 713 &           effective_potential%anharmonics_terms%coefficients(ii)%terms,&
 714 &           effective_potential%anharmonics_terms%coefficients(ii)%name,&
 715 &           check=.false.)
 716          end do
 717 
 718          do ii=1,ncoeff_bound
 719            call polynomial_coeff_init(coeffs_bound(ii)%coefficient,&
 720 &           coeffs_bound(ii)%nterm,&
 721 &           coeffs_all(ncoeff+ii),&
 722 &           coeffs_bound(ii)%terms,&
 723 &           coeffs_bound(ii)%name,&
 724 &           check=.false.)
 725          end do
 726 
 727 !      Copy the fixed coefficients from the model (without bound coeff)
 728          ncoeff = effective_potential%anharmonics_terms%ncoeff
 729          ABI_DATATYPE_ALLOCATE(coeffs_tmp,(ncoeff+ncoeff_bound))
 730          do ii=1,ncoeff
 731            call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
 732 &           effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
 733 &           coeffs_tmp(ii),&
 734 &           effective_potential%anharmonics_terms%coefficients(ii)%terms,&
 735 &           effective_potential%anharmonics_terms%coefficients(ii)%name,&
 736 &           check=.false.)
 737          end do
 738 
 739          ncoeff_max = ncoeff+ncoeff_bound
 740          ABI_ALLOCATE(listcoeff,(ncoeff_max))
 741          listcoeff = 0
 742          do jj=1,ncoeff
 743            listcoeff(jj) = jj
 744          end do
 745 
 746          model_bound = 0
 747          model_ncoeffbound = 0
 748 
 749          do ii=2,inp%bound_maxCoeff
 750 !        Compute the number of possible combination
 751            nmodels = 1
 752            ABI_ALLOCATE(list_bound,(nmodels,ii))
 753            ABI_ALLOCATE(list_tmp,(ii))
 754            list_bound = 0; list_tmp = 0; kk = 0;  jj = 1
 755 
 756 !        Generate the list of possible combinaison 1st count
 757            call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.false.)
 758            nmodels = kk
 759 
 760            write(message, '(5a,I0,a,I0,a)')ch10,'--',ch10,' Try to bound the model ',&
 761 &           'with ', ii,' additional positive terms (',nmodels,') possibilities'
 762            call wrtout(std_out,message,'COLL')
 763 
 764 !        allocate and generate combinaisons
 765            ABI_DEALLOCATE(list_bound)
 766            ABI_DEALLOCATE(list_tmp)
 767            ABI_ALLOCATE(coeff_values,(nmodels,ncoeff+ii))
 768            ABI_ALLOCATE(listcoeff_bound,(nmodels,ncoeff+ii))
 769            ABI_ALLOCATE(list_bound,(nmodels,ii))
 770            ABI_ALLOCATE(list_tmp,(ii))
 771            ABI_ALLOCATE(isPositive,(nmodels))
 772            list_bound = 0;  listcoeff_bound = 0;  list_tmp = 0; isPositive = 0; kk = 0; jj = 1
 773            call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.true.)
 774 !        Generate the models
 775            do jj=1,nmodels
 776              listcoeff_bound(jj,1:ncoeff) = listcoeff(1:ncoeff)
 777              listcoeff_bound(jj,ncoeff+1:ncoeff+ii) = list_bound(jj,:) + ncoeff
 778            end do
 779 
 780 !        Reset the simulation
 781            call effective_potential_setCoeffs(coeffs_all,effective_potential,ncoeff+ncoeff_bound)
 782            call fit_polynomial_coeff_getPositive(effective_potential,hist,coeff_values,&
 783 &           isPositive,listcoeff_bound,ncoeff+ii,ncoeff,nmodels,comm,verbose=.false.)
 784            if(all(isPositive == 0)) then
 785              write(message, '(5a,I0,a)')ch10,'--',ch10,' No possible model ',&
 786 &             'with ', ii,' additional terms found'
 787              call wrtout(std_out,message,'COLL')
 788            else
 789 
 790              do jj=1,nmodels
 791                if(isPositive(jj) == 1 .and. all(abs(coeff_values(jj,:)) < 1.0E5)) then
 792                  write(message, '(2a,I0,a)') ch10,' The model number ',jj,' ['
 793                  do kk=1,ncoeff+ii
 794                    if(kk<ncoeff+ii)then
 795                      write(message, '(a,I0,a)') trim(message),listcoeff_bound(jj,kk),','
 796                    else
 797                      write(message, '(a,I0)') trim(message),listcoeff_bound(jj,kk)
 798                    end if
 799                  end do
 800                  write(message, '(2a)') trim(message),'] is positive'
 801                  call wrtout(std_out,message,'COLL')
 802                  write(message, '(2a,I0,a)') ' Check if the model ',&
 803 &                 'number ', jj,' is bounded...'
 804                  call wrtout(std_out,message,'COLL')
 805 
 806 !             Set the coefficients of the model
 807                  do kk=1,ncoeff+ii
 808                    if(kk<=ncoeff)then
 809 !                 just set the values of the coefficient
 810                      call polynomial_coeff_setCoefficient(coeff_values(jj,kk),coeffs_tmp(kk))
 811                      write(message, '(a,I0,a,ES19.10,2a)') ' Set the value of the coefficient ',kk,&
 812 &                     ' =>',coeff_values(jj,kk),'     ',trim(coeffs_tmp(kk)%name)
 813                      call wrtout(std_out,message,'COLL')
 814 
 815                    else
 816 !                 Set the good coefficient
 817                      icoeff_bound = listcoeff_bound(jj,kk)-ncoeff ! need to remove ncoeff value
 818                      write(message, '(a,I0,a,I0,a,ES19.10,2a)')&
 819 &                     ' Set the value of the coefficient ',kk,' (',icoeff_bound,&
 820 &                     ') =>',coeff_values(jj,kk),&
 821 &                     '     ',trim(coeffs_bound(icoeff_bound)%name)
 822                      call wrtout(std_out,message,'COLL')
 823                      call polynomial_coeff_free(coeffs_tmp(kk))
 824                      call polynomial_coeff_init(coeff_values(jj,kk),&
 825 &                     coeffs_bound(icoeff_bound)%nterm,&
 826 &                     coeffs_tmp(kk),&
 827 &                     coeffs_bound(icoeff_bound)%terms,&
 828 &                     coeffs_bound(icoeff_bound)%name,&
 829 &                     check=.false.)
 830 
 831                    end if
 832                  end do
 833 
 834 !              Reset the simulation and set the coefficients of the model
 835                  call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,&
 836 &                 ncoeff+ii)
 837                  call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),1,0,&
 838 &                 -1,1,comm,verbose=.false.,positive=.false.)
 839                  call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
 840                  dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
 841                  acell(:) = one
 842                  call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
 843 &                 effective_potential%supercell%xcart,xred)
 844                  xred_old = xred
 845                  vel_cell(:,:) = zero
 846                  vel(:,:)      = zero
 847                  fred(:,:)     = zero
 848                  fcart(:,:)    = zero
 849 
 850 !              Run mover
 851                  call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
 852 &                 rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
 853 &                 effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST)
 854 
 855                  if(.not.effective_potential%anharmonics_terms%bounded)then
 856                    write(message, '(2a)' ) ' => The model is not bounded'
 857                  else
 858                    write(message, '(2a)' ) ' => The model is bounded'
 859                  end if
 860                  call wrtout(std_out,message,'COLL')
 861 !             Exit if the model is bounded
 862                  if(effective_potential%anharmonics_terms%bounded) then
 863                    model_bound = jj
 864                    model_ncoeffbound = ii
 865                    exit
 866                  end if
 867                end if
 868              end do
 869            end if
 870 
 871            ABI_DEALLOCATE(list_tmp)
 872            ABI_DEALLOCATE(list_bound)
 873            ABI_DEALLOCATE(isPositive)
 874 
 875 !        Exit if the model is bounded
 876            if(effective_potential%anharmonics_terms%bounded) then
 877 !         Final transfert
 878              write(message, '(3a)' ) ch10,' => The model is now bounded'
 879              call wrtout(ab_out,message,'COLL')
 880              call wrtout(std_out,message,'COLL')
 881              do kk=ncoeff+1,ncoeff+model_ncoeffbound
 882                icoeff_bound = listcoeff_bound(model_bound,kk)-ncoeff ! need to remove ncoeff value
 883                call polynomial_coeff_free(coeffs_tmp(kk))
 884                call polynomial_coeff_init(coeff_values(model_bound,kk),&
 885 &               coeffs_bound(icoeff_bound)%nterm,&
 886 &               coeffs_tmp(kk),&
 887 &               coeffs_bound(icoeff_bound)%terms,&
 888 &               coeffs_bound(icoeff_bound)%name,&
 889 &               check=.false.)
 890              end do
 891              ABI_DEALLOCATE(coeff_values)
 892              ABI_DEALLOCATE(listcoeff_bound)
 893              exit
 894            end if
 895            ABI_DEALLOCATE(coeff_values)
 896 
 897          end do
 898 
 899          do ii=1,ncoeff_bound
 900            call polynomial_coeff_free(coeffs_bound(ii))
 901          end do
 902          if(allocated(coeffs_bound)) ABI_DEALLOCATE(coeffs_bound)
 903 
 904        end if
 905 
 906        if(.not.effective_potential%anharmonics_terms%bounded)then
 907          write(message, '(3a)' ) ch10,' => The model cannot be bounded'
 908          call wrtout(ab_out,message,'COLL')
 909          call wrtout(std_out,message,'COLL')
 910          model_ncoeffbound = 0
 911          model_bound = 0
 912        end if
 913 
 914 !      Fit the final model
 915        call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+model_ncoeffbound),effective_potential,&
 916 &       ncoeff+model_ncoeffbound)
 917 
 918        call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,&
 919 &       -1,1,comm,verbose=.false.)
 920 
 921        write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit bound process: '
 922        call wrtout(ab_out,message,'COLL')
 923        call wrtout(std_out,message,'COLL')
 924 
 925        do ii = 1,ncoeff+model_ncoeffbound
 926          write(message, '(a,I0,a,ES19.10,2a)') " ",ii," =>",&
 927 &         effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
 928 &         " ",trim(effective_potential%anharmonics_terms%coefficients(ii)%name)
 929          call wrtout(ab_out,message,'COLL')
 930          call wrtout(std_out,message,'COLL')
 931        end do
 932 
 933 !      Deallocation
 934        ABI_DEALLOCATE(listcoeff)
 935        do ii=1,ncoeff_max
 936          call polynomial_coeff_free(coeffs_tmp(ii))
 937        end do
 938        if(allocated(coeffs_tmp)) ABI_DEALLOCATE(coeffs_tmp)
 939 
 940        do ii=1,ncoeff_max
 941          call polynomial_coeff_free(coeffs_all(ii))
 942        end do
 943        if(allocated(coeffs_all)) ABI_DEALLOCATE(coeffs_all)
 944 
 945      end if
 946 
 947    else  if (option == -3) then
 948 
 949 !*************************************************************
 950 !   Call the routine for calculation of the energy for specific
 951 !   partern of displacement or strain for the effective
 952 !   Hamiltonian
 953 !*************************************************************
 954 !     write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,&
 955 !&     ' Effective Hamiltonian calculation'
 956 !     call wrtout(ab_out,message,'COLL')
 957 !     call wrtout(std_out,message,'COLL')
 958 
 959 !     acell = one
 960 !     call gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,&
 961 !&                mpi_enreg,npwtot,dtset%occ_orig,pawang,pawrad,pawtab,&
 962 !&                psps,results_gs,dtset%rprimd_orig,scf_history,vel,vel_cell,wvl,xred)
 963 
 964    end if
 965 
 966 !***************************************************************
 967 ! 5   Deallocation of array
 968 !***************************************************************
 969 
 970    ABI_DEALLOCATE(fred)
 971    ABI_DEALLOCATE(fcart)
 972    ABI_DEALLOCATE(indsym)
 973    ABI_DEALLOCATE(rhog)
 974    ABI_DEALLOCATE(rhor)
 975    ABI_DEALLOCATE(vel)
 976    ABI_DEALLOCATE(xred)
 977    ABI_DEALLOCATE(xred_old)
 978    ABI_DEALLOCATE(ab_xfh%xfhist)
 979 
 980    ! if(option == -3)then
 981    !   call args_gs_free(args_gs)
 982    !   call psps_free(psps)
 983    !   do ii = 1,npsp
 984    !     call paw_setup_free(paw_setup(ii))
 985    !   end do
 986    !   ABI_DEALLOCATE(paw_setup)
 987    !   ABI_DEALLOCATE(ipsp2xml)
 988    !   ABI_DEALLOCATE(pspheads)
 989    !   call pawrad_free(pawrad)
 990    !   call pawtab_free(pawtab)
 991    !   ABI_DATATYPE_DEALLOCATE(pawrad)
 992    !   ABI_DATATYPE_DEALLOCATE(pawtab)
 993    !   ABI_DEALLOCATE(npwtot)
 994    ! end if
 995    call dtset_free(dtset)
 996    call destroy_results_gs(results_gs)
 997    call scfcv_destroy(scfcv_args)
 998    call destroy_mpi_enreg(mpi_enreg)
 999 
1000  end if
1001 
1002  write(message, '(a,(80a),a,a)' ) ch10,&
1003 & ('=',ii=1,80),ch10
1004  call wrtout(ab_out,message,'COLL')
1005  call wrtout(std_out,message,'COLL')
1006 
1007 end subroutine mover_effpot