TABLE OF CONTENTS


ABINIT/erlxconv [ Functions ]

[ Top ] [ Functions ]

NAME

  erlxconv

FUNCTION

  FIXME: add description.

INPUTS

  argin(sizein)=description

OUTPUT

  argout(sizeout)=description

PARENTS

      mover

CHILDREN

      wrtout

SOURCE

1146 subroutine erlxconv(hist,iexit,itime,itime_hist,ntime,tolmxde)
1147 
1148 
1149 !This section has been created automatically by the script Abilint (TD).
1150 !Do not modify the following lines by hand.
1151 #undef ABI_FUNC
1152 #define ABI_FUNC 'erlxconv'
1153 !End of the abilint section
1154 
1155  implicit none
1156 
1157 !Arguments ------------------------------------
1158 !scalars
1159  integer,intent(in) :: itime,itime_hist,ntime
1160  integer,intent(inout) :: iexit
1161  real(dp), intent(in) :: tolmxde
1162 !arrays
1163  type(abihist),intent(inout) :: hist
1164 
1165 !Local variables-------------------------------
1166  integer :: ihist,ihist_prev,ihist_prev2
1167  real(dp) :: ediff1,ediff2,maxediff
1168  character(len=500) :: message
1169 ! *************************************************************************
1170 
1171  if (itime_hist<3) then
1172    write(message, '(a,a,a)' ) ch10,&
1173 &   ' erlxconv : minimum 3 Broyd/MD steps to check convergence of energy in relaxations',ch10
1174    call wrtout(std_out,message,'COLL')
1175  else
1176    ihist = hist%ihist
1177    ihist_prev  = abihist_findIndex(hist,-1)
1178    ihist_prev2 = abihist_findIndex(hist,-2)
1179    ediff1 = hist%etot(ihist) - hist%etot(ihist_prev)
1180    ediff2 = hist%etot(ihist) - hist%etot(ihist_prev2)
1181    if ((abs(ediff1)<tolmxde).and.(abs(ediff2)<tolmxde)) then
1182      write(message, '(a,a,i4,a,a,a,a,a,es11.4,a,a)' ) ch10,&
1183 &     ' At Broyd/MD step',itime,', energy is converged : ',ch10,&
1184 &     '  the difference in energy with respect to the two ',ch10,&
1185 &     '  previous steps is < tolmxde=',tolmxde,' ha',ch10
1186      call wrtout(ab_out,message,'COLL')
1187      call wrtout(std_out,message,'COLL')
1188      iexit=1
1189    else
1190      maxediff = max(abs(ediff1),abs(ediff2))
1191      if(iexit==1)then
1192        write(message, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1193 &       ' erlxconv : WARNING -',ch10,&
1194 &       '  ntime=',ntime,' was not enough Broyd/MD steps to converge energy: ',ch10,&
1195 &       '  max difference in energy =',maxediff,' > tolmxde=',tolmxde,' ha',ch10
1196        call wrtout(std_out,message,'COLL')
1197        call wrtout(ab_out,message,'COLL')
1198 
1199        write(std_out,"(8a)")ch10,&
1200 &       "--- !RelaxConvergenceWarning",ch10,&
1201 &       "message: | ",ch10,TRIM(indent(message)),ch10,&
1202 &       "..."
1203      else
1204        write(message, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1205 &       ' erlxconv : at Broyd/MD step',itime,', energy has not converged yet. ',ch10,&
1206 &       '  max difference in energy=',maxediff,' > tolmxde=',tolmxde,' ha',ch10
1207        call wrtout(std_out,message,'COLL')
1208      end if
1209    end if
1210  end if
1211 
1212 end subroutine erlxconv

ABINIT/fconv [ Functions ]

[ Top ] [ Functions ]

NAME

 fconv

FUNCTION

 Check maximal absolute value of force (hartree/bohr) against
 input tolerance; if below tolerance, return iexit=1.
 Takes into account the fact that the Broyden (or moldyn) step
 might be the last one (last itime), to print eventually modified message.
 Stresses are also included in the check, provided that optcell/=0.
 If optcell=1, takes only the trace into account
    optcell=2, takes all components into account
    optcell=3, takes traceless stress into account
    optcell=4, takes sigma(1 1) into account
    optcell=5, takes sigma(2 2) into account
    optcell=6, takes sigma(3 3) into account
    optcell=7, takes sigma(2,2),(2,3) and (3 3) into account
    optcell=8, takes sigma(1,1),(1,3) and (3 3) into account
    optcell=9, takes sigma(1,1),(1,2) and (2 2) into account
 In the case of stresses, target the tensor strtarget, and
 take into account the factor strfact

INPUTS

  fcart(3,natom)= forces on atoms in hartree/bohr in cartesian coordinates
  iatfix(3,natom)=1 for frozen atom, 0 for unfrozen
  itime=current number of Broyden/Moldyn iterations
  natom=number of atoms in unit cell
  ntime=maximum number of Broyden/Moldyn iterations allowed
  optcell=option for taking stresses into account (see above)
  strfact=factor that multiplies the stresses when they are compared to forces.
  strtarget(6)=components of the target stress tensor (hartree/bohr^3)
  strten(6)=components of the stress tensor (hartree/bohr^3)
  tolmxf=tolerance on maximal absolute value of components of forces

OUTPUT

  writes to unit std_out and to ab_out, and returns

SIDE EFFECTS

 Input/Output
  at input  : iexit=  0 if not the last itime,  1 if the last itime
  at output : iexit=  0 if not below tolerance, 1 if below tolerance

PARENTS

      mover

CHILDREN

      wrtout

SOURCE

1017 subroutine fconv(fcart,iatfix,iexit,itime,natom,ntime,optcell,strfact,strtarget,strten,tolmxf)
1018 
1019 
1020 !This section has been created automatically by the script Abilint (TD).
1021 !Do not modify the following lines by hand.
1022 #undef ABI_FUNC
1023 #define ABI_FUNC 'fconv'
1024 !End of the abilint section
1025 
1026  implicit none
1027 
1028 !Arguments ------------------------------------
1029 !scalars
1030  integer,intent(in) :: itime,natom,ntime,optcell
1031  integer,intent(inout) :: iexit
1032  real(dp),intent(in) :: strfact,tolmxf
1033 !arrays
1034  integer,intent(in) :: iatfix(3,natom)
1035  real(dp),intent(in) :: fcart(3,natom),strtarget(6),strten(6)
1036 
1037 !Local variables-------------------------------
1038 !scalars
1039  integer :: iatom,idir,istr
1040  real(dp) :: fmax,strdiag
1041  character(len=500) :: message
1042 !arrays
1043  real(dp) :: dstr(6)
1044 
1045 ! *************************************************************************
1046 
1047 !Compute maximal component of forces, EXCLUDING any fixed components
1048  fmax=zero
1049  do iatom=1,natom
1050    do idir=1,3
1051      if (iatfix(idir,iatom) /= 1) then
1052        if( abs(fcart(idir,iatom)) >= fmax ) fmax=abs(fcart(idir,iatom))
1053      end if
1054    end do
1055  end do
1056 
1057  dstr(:)=strten(:)-strtarget(:)
1058 
1059 !Eventually take into account the stress
1060  if(optcell==1)then
1061    strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp
1062    if(abs(strdiag)*strfact >= fmax ) fmax=abs(strdiag)*strfact
1063  else if(optcell==2)then
1064    do istr=1,6
1065      if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact
1066    end do
1067  else if(optcell==3)then
1068 !  Must take away the trace from diagonal elements
1069    strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp
1070    do istr=1,3
1071      if(abs(dstr(istr)-strdiag)*strfact >= fmax ) fmax=abs(dstr(istr)-strdiag)*strfact
1072    end do
1073    do istr=4,6
1074      if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact
1075    end do
1076  else if(optcell==4 .or. optcell==5 .or. optcell==6)then
1077    if(abs(dstr(optcell-3))*strfact >= fmax ) fmax=abs(dstr(optcell-3))*strfact
1078  else if(optcell==7)then
1079    if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact
1080    if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact
1081    if(abs(dstr(4))*strfact >= fmax ) fmax=abs(dstr(4))*strfact
1082  else if(optcell==8)then
1083    if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact
1084    if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact
1085    if(abs(dstr(5))*strfact >= fmax ) fmax=abs(dstr(5))*strfact
1086  else if(optcell==9)then
1087    if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact
1088    if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact
1089    if(abs(dstr(6))*strfact >= fmax ) fmax=abs(dstr(6))*strfact
1090  end if
1091 
1092  if (fmax<tolmxf) then
1093    write(message, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1094 &   ' At Broyd/MD step',itime,', gradients are converged : ',ch10,&
1095 &   '  max grad (force/stress) =',fmax,' < tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1096    call wrtout(ab_out,message,'COLL')
1097    call wrtout(std_out,message,'COLL')
1098    iexit=1
1099  else
1100    if(iexit==1)then
1101      write(message, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1102 &     ' fconv : WARNING -',ch10,&
1103 &     '  ntime=',ntime,' was not enough Broyd/MD steps to converge gradients: ',ch10,&
1104 &     '  max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1105      call wrtout(std_out,message,'COLL')
1106      call wrtout(ab_out,message,'COLL')
1107 
1108      write(std_out,"(8a)")ch10,&
1109 &     "--- !RelaxConvergenceWarning",ch10,&
1110 &     "message: | ",ch10,TRIM(indent(message)),ch10,&
1111 &     "..."
1112 
1113    else
1114      write(message, '(a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) &
1115 &     ' fconv : at Broyd/MD step',itime,', gradients have not converged yet. ',ch10,&
1116 &     '  max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1117      call wrtout(std_out,message,'COLL')
1118    end if
1119    iexit=0
1120  end if
1121 
1122 end subroutine fconv

ABINIT/gettag [ Functions ]

[ Top ] [ Functions ]

NAME

 gettag

FUNCTION

 Set the tag associated to each atom,

INPUTS

 prtallatoms = Logical for PRTint ALL ATOMS
 atlist      = ATom LIST
 index       = index for each atom
 natom       = Number of ATOMs

OUTPUT

  tag = The string to put for each atom

PARENTS

      prtxfase

CHILDREN

      gettag,wrtout

SOURCE

1526 subroutine gettag(atlist,index,natom,prtallatoms,tag)
1527 
1528 
1529 !This section has been created automatically by the script Abilint (TD).
1530 !Do not modify the following lines by hand.
1531 #undef ABI_FUNC
1532 #define ABI_FUNC 'gettag'
1533 !End of the abilint section
1534 
1535 implicit none
1536 
1537 !Arguments ------------------------------------
1538 !scalars
1539   logical,intent(in) :: prtallatoms
1540   integer,intent(in) :: natom
1541   logical,intent(in) :: atlist(natom)
1542   integer,intent(in) :: index
1543   character(len=7),intent(out)   :: tag
1544 
1545 !Local variables -------------------------
1546 
1547 ! *********************************************************************
1548 !The numbering will be from (1) to (9999)
1549 
1550  if (prtallatoms)then
1551    tag=''
1552  elseif (atlist(index)) then
1553    if (natom<10) then
1554      write(tag, '(a,I1.1,a)') ' (',index,')'
1555    elseif (natom<100) then
1556      write(tag, '(a,I2.2,a)') ' (',index,')'
1557    elseif (natom<1000) then
1558      write(tag, '(a,I3.3,a)') ' (',index,')'
1559    elseif (natom<10000) then
1560      write(tag, '(a,I4.4,a)') ' (',index,')'
1561    end if
1562  end if
1563 
1564  end subroutine gettag

ABINIT/m_mover [ Modules ]

[ Top ] [ Modules ]

NAME

  m_mover

FUNCTION

 Move ion or change acell according to forces and stresses

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, SE)
  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
28 
29  use defs_basis
30  use defs_abitypes
31  use m_errors
32  use m_abimover
33  use m_abihist
34  use m_xmpi
35  use m_nctk
36 #ifdef HAVE_NETCDF
37  use netcdf
38 #endif
39 #if defined HAVE_LOTF
40  use lotfpath
41  use m_pred_lotf
42 #endif
43 
44  use m_fstrings,           only : strcat, sjoin, indent
45  use m_symtk,              only : matr3inv, symmetrize_xred
46  use m_geometry,           only : fcart2fred, chkdilatmx, xred2xcart
47  use m_crystal,            only : crystal_init, crystal_free, crystal_t
48  use m_crystal_io,         only : crystal_ncwrite_path
49  use m_time,               only : abi_wtime, sec2str
50  use m_exit,               only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
51  use m_electronpositron,   only : electronpositron_type
52  use m_scfcv,              only : scfcv_t, scfcv_run
53  use m_effective_potential,only : effective_potential_type,effective_potential_evaluate
54  use m_dtfil,              only : dtfil_init_time, status
55  use m_initylmg,           only : initylmg
56  use m_xfpack,             only : xfh_update
57  use m_precpred_1geo,      only : precpred_1geo
58  use m_pred_delocint,      only : pred_delocint
59  use m_pred_bfgs,          only : pred_bfgs, pred_lbfgs
60  use m_pred_fire,          only : pred_fire
61  use m_pred_isokinetic,    only : pred_isokinetic
62  use m_pred_diisrelax,     only : pred_diisrelax
63  use m_pred_nose,          only : pred_nose
64  use m_pred_srkhna14,      only : pred_srkna14
65  use m_pred_isothermal,    only : pred_isothermal
66  use m_pred_verlet,        only : pred_verlet
67  use m_pred_velverlet,     only : pred_velverlet
68  use m_pred_moldyn,        only : pred_moldyn
69  use m_pred_langevin,      only : pred_langevin
70  use m_pred_steepdesc,     only : pred_steepdesc
71  use m_pred_simple,        only : pred_simple, prec_simple
72  use m_pred_hmc,           only : pred_hmc
73  use m_generate_training_set, only : generate_training_set
74  use m_wvl_wfsinp, only : wvl_wfsinp_reformat
75  use m_wvl_rho,      only : wvl_mkrho
76 
77  implicit none
78 
79  private

ABINIT/mover [ Functions ]

[ Top ] [ Functions ]

NAME

 mover

FUNCTION

 Move ion or change acell acording to forces and stresses

INPUTS

  amu_curr(ntypat)=mass of each atom for the current image
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs
   | mkmem =number of k points treated by this node
   |  angular momentum for nonlocal pseudopotential
   | mpw=maximum dimensioned size of npw.
   | natom=number of atoms in unit cell
   |  except on first call (hartree/bohr); updated on output
   | nfft=(effective) number of FFT grid points (for this processor)
   |      for the "coarse" grid (see NOTES below)
   | nkpt=number of k points.
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in space group
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mpi_enreg=informations about MPI parallelization
  nfftf=(effective) number of FFT grid points (for this processor)
       for the "fine" grid (see NOTES below)
  npwarr(nkpt)=number of planewaves in basis and boundary at this k point.
  nattyp(ntypat)= # atoms of each type.
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
   | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  rprimd(3,3)=dimensional primitive translations (bohr)

OUTPUT

  results_gs <type(results_gs_type)>=results (energy and its components,
   forces and its components, the stress tensor) of a ground-state computation
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  resid(mband*nkpt*nsppol)=residuals for each band over all k points.

SIDE EFFECTS

 Rest of i/o is related to lda
  acell(3)=length scales of primitive translations (bohr)
  cg(2,mcg)=array for planewave coefficients of wavefunctions.
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  initialized= if 0 the initialisation of the gstate run is not yet finished
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
  occ(mband*nkpt*nsppol=occupation number for each band (usually 2) at each k point.
  rhog(2,nfftf)=array for Fourier transform of electron density
  rhor(nfftf,nspden)=array for electron density in electrons/bohr**3.
  scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
  taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
  vel(3,natom)=old value of velocity; updated on output
  vel_cell(3,3)=old value of cell parameters velocity; updated on output
  xred(3,natom)=reduced dimensionless atomic coordinates; updated on output
  xred_old(3,natom)=work space for old xred
  eff_pot<type(effective_potential_type)> = optional,effective_potential datatype
  verbose = optional, default is true, flag to disable the verbose mode
  write_HIST = optional, default is true, flag to disble the write of the HIST file

NOTES

 This subroutine uses the arguments natom, xred, vel, amu_curr,
 vis, and dtion (the last two contained in dtset) to make
 molecular dynamics updates.  The rest of the lengthy
 argument list supports the underlying lda computation
 of forces, returned from subroutine scfcv

 USE OF FFT GRIDS:
 =================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut)
      for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ...
      are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg)
      for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ...
      Total density, potentials, ...
      are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used.
      It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf)
      are set equal to (nfft,ngfft,mgfft) in that case.

PARENTS

      gstate,mover_effpot

CHILDREN

      abiforstr_fin,abiforstr_ini,abihist_bcast,abihist_compare_and_copy
      abihist_free,abihist_init,abimover_fin,abimover_ini
      chkdilatmx,crystal_free,crystal_init,dtfil_init_time
      effective_potential_evaluate,erlxconv,fcart2fred,fconv,hist2var
      initylmg,matr3inv,mttk_fin,mttk_ini,prec_simple,pred_bfgs,pred_delocint
      pred_diisrelax,pred_hmc,pred_isokinetic,pred_isothermal,pred_langevin
      pred_lbfgs,pred_lotf,pred_moldyn,pred_nose,pred_simple,pred_srkna14
      pred_steepdesc,pred_velverlet,pred_verlet,prtxfase,read_md_hist
      scfcv_run,status,symmetrize_xred,var2hist,vel2hist,write_md_hist
      wrt_moldyn_netcdf,wrtout,wvl_mkrho,wvl_wfsinp_reformat,xfh_update
      xmpi_barrier,xmpi_isum,xmpi_wait

SOURCE

199 subroutine mover(scfcv_args,ab_xfh,acell,amu_curr,dtfil,&
200 & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,&
201 & effective_potential,filename_ddb,verbose,writeHIST)
202 
203 
204 !This section has been created automatically by the script Abilint (TD).
205 !Do not modify the following lines by hand.
206 #undef ABI_FUNC
207 #define ABI_FUNC 'mover'
208 !End of the abilint section
209 
210 implicit none
211 
212 !Arguments ------------------------------------
213 !scalars
214 type(scfcv_t),intent(inout) :: scfcv_args
215 type(datafiles_type),intent(inout),target :: dtfil
216 type(electronpositron_type),pointer :: electronpositron
217 type(ab_xfh_type),intent(inout) :: ab_xfh
218 type(effective_potential_type),optional,intent(inout) :: effective_potential
219 logical,optional,intent(in) :: verbose
220 logical,optional,intent(in) :: writeHIST
221 character(len=fnlen),optional,intent(in) :: filename_ddb
222 !arrays
223 real(dp),intent(inout) :: acell(3)
224 real(dp), intent(in),target :: amu_curr(:) !(scfcv%dtset%ntypat)
225 real(dp), pointer :: rhog(:,:),rhor(:,:)
226 real(dp), intent(inout) :: xred(3,scfcv_args%dtset%natom),xred_old(3,scfcv_args%dtset%natom)
227 real(dp), intent(inout) :: vel(3,scfcv_args%dtset%natom),vel_cell(3,3),rprimd(3,3)
228 
229 !Local variables-------------------------------
230 !scalars
231 integer,parameter :: level=102,master=0
232 type(abihist) :: hist,hist_prev
233 type(abimover) :: ab_mover
234 type(abimover_specs) :: specs
235 type(abiforstr) :: preconforstr ! Preconditioned forces and stress
236 type(delocint) :: deloc
237 type(mttk_type) :: mttk_vars
238 integer :: irshift,itime,icycle,itime_hist,iexit=0,ifirst,ihist_prev,ihist_prev2,timelimit_exit,ncycle,nhisttot,kk,jj,me
239 integer :: nloop,nshell,ntime,option,comm
240 integer :: nerr_dilatmx,my_quit,ierr,quitsum_request
241 integer ABI_ASYNC :: quitsum_async
242 character(len=500) :: message
243 character(len=500) :: dilatmx_errmsg
244 character(len=8) :: stat4xml
245 character(len=35) :: fmt
246 character(len=fnlen) :: filename,fname_ddb
247 real(dp) :: favg
248 logical :: DEBUG=.FALSE., need_verbose=.TRUE.,need_writeHIST=.TRUE.
249 logical :: need_scfcv_cycle = .TRUE.
250 logical :: change,useprtxfase
251 logical :: skipcycle
252 integer :: minIndex,ii,similar,conv_retcode
253 integer :: iapp
254 real(dp) :: minE,wtime_step,now,prev
255 type(crystal_t) :: crystal
256 logical :: file_exists
257 !arrays
258 real(dp) :: gprimd(3,3),rprim(3,3),rprimd_prev(3,3)
259 real(dp),allocatable :: fred_corrected(:,:),xred_prev(:,:)
260 ! ***************************************************************
261  need_verbose=.TRUE.
262  if(present(verbose)) need_verbose = verbose
263 
264  need_writeHIST=.TRUE.
265  if(present(writeHIST)) need_writeHIST = writeHIST
266 
267  call status(0,dtfil%filstat,iexit,level,'init          ')
268 
269  ! enable time limit handler if not done in callers.
270  if (enable_timelimit_in(ABI_FUNC) == ABI_FUNC) then
271    if (need_verbose) then
272      write(std_out,*)"Enabling timelimit check in function: ",trim(ABI_FUNC)," with timelimit: ",trim(sec2str(get_timelimit()))
273    end if
274  end if
275 
276 !Table of contents
277 !(=>) Refers to an important call (scfcv,pred_*)
278 !
279 !01. Initialization of indexes and allocations of arrays
280 !02. Particularities of each predictor
281 !03. Set the number of iterations ntime
282 !04. Try to read history of previous calculations
283 !05. Allocate the hist structure
284 !06. First output before any itime or icycle
285 !07. Fill the history of the first SCFCV
286 !08. Loop for itime (From 1 to ntime)
287 !09. Loop for icycle (From 1 to ncycle)
288 !10. Output for each icycle (and itime)
289 !11. Symmetrize atomic coordinates over space group elements
290 !12. => Call to SCFCV routine and fill history with forces
291 !13. Write the history into the _HIST file
292 !14. Output after SCFCV
293 !15. => Test Convergence of forces and stresses
294 !16. => Precondition forces, stress and energy
295 !17. => Call to each predictor
296 !18. Use the history  to extract the new values
297 !19. End loop icycle
298 !20. End loop itime
299 !21. Set the final values of xred
300 !22. XML Output at the end
301 !23. Deallocate hist and ab_mover datatypes
302 !
303  call abimover_ini(ab_mover,amu_curr,dtfil,scfcv_args%dtset,specs)
304 
305  if(ab_mover%ionmov==10 .or. ab_mover%ionmov==11)then
306    call delocint_ini(deloc)
307  end if
308 
309  if (ab_mover%ionmov==13)then
310    call mttk_ini(mttk_vars,ab_mover%nnos)
311  end if
312 
313 !###########################################################
314 !### 03. Set the number of iterations ntime
315 !###     By default ntime==1 but if the user enter a lower
316 !###     value mover will execute at least one iteration
317 
318  if (scfcv_args%dtset%ntime<1)then
319    ntime=1
320  else
321    ntime=scfcv_args%dtset%ntime
322  end if
323 
324 !###########################################################
325 !### 04. Try to read history of previous calculations
326 !###     It requires access to the NetCDF library
327 
328 !Init MPI data
329  comm=scfcv_args%mpi_enreg%comm_cell
330  me=xmpi_comm_rank(comm)
331 
332 #if defined HAVE_NETCDF
333  filename=trim(ab_mover%filnam_ds(4))//'_HIST.nc'
334 
335  if (ab_mover%restartxf<0)then
336 !  Read history from file (and broadcast if MPI)
337    if (me==master) then
338      call read_md_hist(filename,hist_prev,specs%isVused,specs%isARused,ab_mover%restartxf==-3)
339    end if
340    call abihist_bcast(hist_prev,master,comm)
341 
342 !  If restartxf specifies to reconstruct the history
343    if (hist_prev%mxhist>0.and.ab_mover%restartxf==-1)then
344      ntime=ntime+hist_prev%mxhist
345    end if
346 
347 !  If restartxf specifies to start from the lowest energy
348    if (hist_prev%mxhist>0.and.ab_mover%restartxf==-2)then
349      minE=hist_prev%etot(1)
350      minIndex=1
351      do ii=1,hist_prev%mxhist
352        if(need_verbose) write(std_out,*) 'Iteration:',ii,' Total Energy:',hist_prev%etot(ii)
353        if (minE>hist_prev%etot(ii))then
354          minE=hist_prev%etot(ii)
355          minIndex=ii
356        end if
357      end do
358      if(need_verbose)write(std_out,*) 'The lowest energy occurs at iteration:',minIndex,'etotal=',minE
359      acell(:)   =hist_prev%acell(:,minIndex)
360      rprimd(:,:)=hist_prev%rprimd(:,:,minIndex)
361      xred(:,:)  =hist_prev%xred(:,:,minIndex)
362      call abihist_free(hist_prev)
363    end if
364 !  If restarxf specifies to start to the last iteration
365    if (hist_prev%mxhist>0.and.ab_mover%restartxf==-3)then
366      acell(:)   =hist_prev%acell(:,hist_prev%mxhist)
367      rprimd(:,:)=hist_prev%rprimd(:,:,hist_prev%mxhist)
368      xred(:,:)  =hist_prev%xred(:,:,hist_prev%mxhist)
369      call abihist_free(hist_prev)
370    end if
371 
372  end if !if (ab_mover%restartxf<=0)
373 
374 #endif
375 
376 !###########################################################
377 !### 05. Allocate the hist structure
378 
379  iexit=0; timelimit_exit=0
380  ncycle=specs%ncycle
381 
382  if(ab_mover%ionmov==25.and.scfcv_args%dtset%hmctt>=0) then
383    ncycle=scfcv_args%dtset%hmctt
384    if(scfcv_args%dtset%hmcsst>0.and.ab_mover%optcell/=0) then
385      ncycle=ncycle+scfcv_args%dtset%hmcsst
386    endif
387  endif
388 
389  nhisttot=ncycle*ntime;if (scfcv_args%dtset%nctime>0) nhisttot=nhisttot+1
390 
391 !AM_2017 New version of the hist, we just store the needed history step not all of them...
392  if(specs%nhist/=-1) nhisttot = specs%nhist ! We don't need to store all the history
393  call abihist_init(hist,ab_mover%natom,nhisttot,specs%isVused,specs%isARused)
394  call abiforstr_ini(preconforstr,ab_mover%natom)
395 
396 !###########################################################
397 !### 06. First output before any itime or icycle
398 
399 !If effective potential is present,
400 !  forces will be compute with it
401  if (present(effective_potential)) then
402    need_scfcv_cycle = .FALSE.
403    if(need_verbose)then
404      write(message,'(2a,i2,5a,80a)')&
405 &     ch10,'=== [ionmov=',ab_mover%ionmov,'] ',trim(specs%method),' with effective potential',&
406 &     ch10,('=',kk=1,80)
407      call wrtout(ab_out,message,'COLL')
408      call wrtout(std_out,message,'COLL')
409    end if
410  else
411    if(need_verbose)then
412      write(message,'(a,a,i2,a,a,a,80a)')&
413 &     ch10,'=== [ionmov=',ab_mover%ionmov,'] ',specs%method,&
414 &     ch10,('=',kk=1,80)
415      call wrtout(ab_out,message,'COLL')
416      call wrtout(std_out,message,'COLL')
417    end if
418  end if
419 
420 !Format for printing on each cycle
421  write(fmt,'(a6,i2,a4,i2,a4,i2,a4,i2,a9)')&
422 & '(a,a,i',int(log10(real(ntime))+1),&
423 & ',a,i',int(log10(real(ntime))+1),&
424 & ',a,i',int(log10(real(ncycle))+1),&
425 & ',a,i',int(log10(real(ncycle))+1),&
426 & ',a,a,80a)'
427 
428 !###########################################################
429 !### 07. Fill the history of the first SCFCV
430 
431  if (ab_mover%ionmov==26)then
432 !Tdep call need to merge with adewandre branch
433  else if (ab_mover%ionmov==27)then
434    if(present(filename_ddb))then
435      fname_ddb = trim(filename_ddb)
436    else
437      fname_ddb = trim(ab_mover%filnam_ds(3))//'_DDB'
438    end if
439    INQUIRE(FILE=filename, EXIST=file_exists)
440 
441    call generate_training_set(acell,ab_mover%ph_freez_disp_addStrain==1,ab_mover%ph_freez_disp_ampl,&
442 &                             fname_ddb,hist,ab_mover%natom,ab_mover%ph_freez_disp_nampl,ntime,&
443 &                             ab_mover%ph_ngqpt,ab_mover%ph_nqshift,ab_mover%ph_freez_disp_option,&
444 &                             ab_mover%ph_qshift,scfcv_args%dtset%supercell_latt,&
445 &                             rprimd,ab_mover%mdtemp(2),xred,comm,DEBUG)
446 
447 
448    !Fill history with the values of xred, acell and rprimd of the first configuration
449    acell(:)   =hist%acell(:,1)
450    rprimd(:,:)=hist%rprimd(:,:,1)
451    xred(:,:)  =hist%xred(:,:,1)
452 
453  else
454 
455    !Fill history with the values of xred, acell and rprimd
456    call var2hist(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
457 
458    !Fill velocities and ionic kinetic energy
459    call vel2hist(ab_mover%amass,hist,vel,vel_cell)
460    hist%time(hist%ihist)=zero
461 
462  end if
463 !Decide if prtxfase will be called
464  useprtxfase=.FALSE.
465  do ii=1,ab_mover%natom
466    if (ab_mover%prtatlist(ii)/=0)then
467      useprtxfase=.TRUE.
468      exit
469    end if
470  end do
471 
472 !At beginning no error
473  nerr_dilatmx = 0
474 
475  ABI_ALLOCATE(xred_prev,(3,scfcv_args%dtset%natom))
476 
477 !###########################################################
478 !### 08. Loop for itime (From 1 to ntime)
479  quitsum_request = xmpi_request_null
480 
481  do itime=1,ntime
482 
483    ! Handle time limit condition.
484    if (itime == 1) prev = abi_wtime()
485    if (itime  > 1) then
486      now = abi_wtime()
487      wtime_step = now - prev
488      prev = now
489      write(message,*)sjoin("mover: previous time step took ",sec2str(wtime_step))
490      if(need_verbose)call wrtout(std_out, message, "COLL")
491      if (have_timelimit_in(ABI_FUNC)) then
492        if (itime > 2) then
493          call xmpi_wait(quitsum_request,ierr)
494          if (quitsum_async > 0) then
495            write(message,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),&
496 &           ". Will exit itime loop in mover."
497            if(need_verbose)MSG_COMMENT(message)
498            if(need_verbose)call wrtout(ab_out, message, "COLL")
499            timelimit_exit = 1
500            exit
501          end if
502        end if
503 
504        my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1
505        call xmpi_isum(my_quit,quitsum_async,comm,quitsum_request,ierr)
506      end if
507    end if
508 
509    skipcycle=.FALSE.
510 #if defined HAVE_LOTF
511    if(ab_mover%ionmov==23 .and. .not. lotf_extrapolation(itime)) skipcycle=.True.
512 #endif
513 
514 !  ###########################################################
515 !  ### 09. Loop for icycle (From 1 to ncycle)
516    do icycle=1,ncycle
517 
518      itime_hist = (itime-1)*ncycle + icycle ! Store the time step in of the history
519 
520 !    ###########################################################
521 !    ### 10. Output for each icycle (and itime)
522      if(need_verbose)then
523        write(message,fmt)&
524 &       ch10,'--- Iteration: (',itime,'/',ntime,') Internal Cycle: (',icycle,'/',ncycle,')',ch10,('-',kk=1,80)
525        call wrtout(ab_out,message,'COLL')
526        call wrtout(std_out,message,'COLL')
527      end if
528      if (useprtxfase) then
529        call prtxfase(ab_mover,hist,itime_hist,std_out,mover_BEFORE)
530      end if
531 
532      xred_prev(:,:)=xred(:,:)
533      rprimd_prev(:,:)=rprimd(:,:)
534 
535 !    ###########################################################
536 !    ### 11. Symmetrize atomic coordinates over space group elements
537 
538      call symmetrize_xred(scfcv_args%indsym,ab_mover%natom,&
539 &     scfcv_args%dtset%nsym,scfcv_args%dtset%symrel,scfcv_args%dtset%tnons,xred)
540 
541      change=any(xred(:,:)/=xred_prev(:,:))
542      if (change)then
543        hist%xred(:,:,hist%ihist)=xred(:,:)
544        if(need_verbose) then
545          write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED'
546          write(std_out,*) 'DIFFERENCES:'
547          do kk=1,ab_mover%natom
548            write(std_out,*) xred(:,kk)-xred_prev(:,kk)
549          end do
550        end if
551        xred_prev(:,:)=xred(:,:)
552      end if
553 
554 !    ###########################################################
555 !    ### 12. => Call to SCFCV routine and fill history with forces
556      if (need_verbose) then
557        if (need_scfcv_cycle) then
558          write(message,'(a,3a,33a,44a)')&
559 &         ch10,('-',kk=1,3),&
560 &         'SELF-CONSISTENT-FIELD CONVERGENCE',('-',kk=1,44)
561        else
562          write(message,'(a,3a,33a,44a)')&
563 &         ch10,('-',kk=1,3),&
564 &         'EFFECTIVE POTENTIAL CALCULATION',('-',kk=1,44)
565        end if
566        call wrtout(ab_out,message,'COLL')
567        call wrtout(std_out,message,'COLL')
568      end if
569 
570      if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then
571 
572        call abihist_compare_and_copy(hist_prev,hist,ab_mover%natom,similar,tol8,specs%nhist==nhisttot)
573        hist_prev%ihist=hist_prev%ihist+1
574 
575      else
576        scfcv_args%ndtpawuj=0
577        iapp=itime
578        if(icycle>1.and.icycle/=ncycle) iapp=-1
579        if(itime==1 .and. icycle/=ncycle ) iapp=-icycle-1
580        if (ab_mover%ionmov==14.and.(icycle<ncycle)) iapp=-1
581 
582 #if defined HAVE_LOTF
583        if (ab_mover%ionmov/=23 .or.(lotf_extrapolation(itime).and.(icycle/=1.or.itime==1)))then
584 #endif
585          !call scfcv_new2(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
586 
587          !WVL - reformat the wavefunctions in the case of xred != xred_old
588          if (scfcv_args%dtset%usewvl == 1 .and. maxval(xred_old - xred) > zero) then
589           !Before running scfcv, on non-first geometry step iterations,
590           ! we need to reformat the wavefunctions, taking into acount the new
591           ! coordinates. We prepare to change rhog (to be removed) and rhor.
592            ABI_DEALLOCATE(rhog)
593            ABI_DEALLOCATE(rhor)
594            call wvl_wfsinp_reformat(scfcv_args%dtset, scfcv_args%mpi_enreg,&
595 &           scfcv_args%psps, rprimd, scfcv_args%wvl, xred, xred_old)
596            scfcv_args%nfftf = scfcv_args%dtset%nfft
597            ABI_ALLOCATE(rhog,(2, scfcv_args%dtset%nfft))
598            ABI_ALLOCATE(rhor,(2, scfcv_args%dtset%nfft))
599            call wvl_mkrho(scfcv_args%dtset, scfcv_args%irrzon, scfcv_args%mpi_enreg,&
600 &           scfcv_args%phnons, rhor,scfcv_args%wvl%wfs,scfcv_args%wvl%den)
601          end if
602 
603 !        MAIN CALL TO SELF-CONSISTENT FIELD ROUTINE
604          if (need_scfcv_cycle) then
605 
606            call dtfil_init_time(dtfil,iapp)
607            call scfcv_run(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
608            if (conv_retcode == -1) then
609                message = "Scf cycle returned conv_retcode == -1 (timelimit is approaching), this should not happen inside mover"
610                MSG_WARNING(message)
611            end if
612 
613          else
614 !          For monte carlo don't need to recompute energy here
615 !          (done in pred_montecarlo)
616            call effective_potential_evaluate( &
617 &           effective_potential,scfcv_args%results_gs%etotal,&
618 &           scfcv_args%results_gs%fcart,scfcv_args%results_gs%fred,&
619 &           scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=need_verbose)
620 
621 !          Check if the simulation did not diverge...
622            if(itime > 3 .and.ABS(scfcv_args%results_gs%etotal - hist%etot(1)) > 1E2)then
623 !            We set to false the flag corresponding to the bound
624              effective_potential%anharmonics_terms%bounded = .FALSE.
625              if(need_verbose.and.me==master)then
626                message = "The simulation is diverging, please check your effective potential"
627                MSG_WARNING(message)
628              end if
629 !            Set the flag to finish the simulation
630              iexit=1
631              stat4xml="Failed"
632            else
633 !            We set to true the flag corresponding to the bound
634              effective_potential%anharmonics_terms%bounded = .TRUE.
635            end if
636          end if
637 #if defined HAVE_LOTF
638        end if
639 #endif
640 !      ANOMALOUS SITUATION
641 !      This is the only case where rprimd could change inside scfcv
642 !      It generates an weird condition, we start with a certain
643 !      value for rprimd before scfcv and after we finish with
644 !      a different value.
645 !      Notice that normally scfcv should not change rprimd
646 !      And even worse if optcell==0
647 !      The solution here is to recompute acell and store these value
648 !      in the present record even if initially it was not exactly
649 !      the value entering in scfcv
650 !      One test case with these condition is bigdft/t10
651        if (any(rprimd(:,:)/=rprimd_prev(:,:))) then
652          hist%acell(:,hist%ihist)=acell(:)
653          hist%rprimd(:,:,hist%ihist)=rprimd(:,:)
654        end if
655 
656 !      ANOMALOUS SITUATIONS
657 !      * In ionmov 4 & 5 xred could change inside SCFCV
658 !      So we need to take the values from the output
659 !
660 !      * Inside scfcv_core.F90 there is a call to symmetrize_xred.F90
661 !      for the first SCF cycle symmetrize_xred could change xred
662        if (ab_mover%ionmov<10)then
663          change=any(xred(:,:)/=xred_prev(:,:))
664          if (change)then
665            hist%xred(:,:,hist%ihist)=xred(:,:)
666            if(need_verbose)then
667              write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED AFTER SCFCV'
668              write(std_out,*) 'DIFFERENCES:'
669              do kk=1,ab_mover%natom
670                write(std_out,*) xred(:,kk)-xred_prev(:,kk)
671              end do
672            end if
673          end if !if (change)
674        end if
675 
676 !      Fill velocities and ionic kinetic energy
677        call vel2hist(ab_mover%amass,hist,vel,vel_cell)
678 
679        hist%fcart(:,:,hist%ihist)=scfcv_args%results_gs%fcart(:,:)
680        hist%strten(:,hist%ihist) =scfcv_args%results_gs%strten(:)
681        hist%etot(hist%ihist)     =scfcv_args%results_gs%etotal
682        hist%entropy(hist%ihist)  =scfcv_args%results_gs%energies%entropy
683        hist%time(hist%ihist)     =real(itime,kind=dp)
684 
685 !      !######################################################################
686      end if ! if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then
687 
688 !    Store trajectory in xfh file
689      if((ab_xfh%nxfh==0.or.itime/=1)) then
690        ABI_ALLOCATE(fred_corrected,(3,scfcv_args%dtset%natom))
691        call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom)
692 !      Get rid of mean force on whole unit cell,
693 !       but only if no generalized constraints are in effect
694        if (ab_mover%nconeq==0)then
695          do ii=1,3
696            if (ii/=3.or.ab_mover%jellslab==0) then
697              favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
698              fred_corrected(ii,:)=fred_corrected(ii,:)-favg
699            end if
700          end do
701        end if
702        if (ncycle<10.and.ab_mover%restartxf>=0) then
703          do ii=1,3
704            rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
705          end do
706 
707 !        The size of ab_xfh%xfhist is to big for very large supercell.
708 !        Call it only for specific ionmov
709          if(any((/2,3,10,11,22/)==ab_mover%ionmov)) then
710            call xfh_update(ab_xfh,acell,fred_corrected,ab_mover%natom,rprim,&
711 &           hist%strten(:,hist%ihist),xred)
712          end if
713        end if
714        ABI_DEALLOCATE(fred_corrected)
715      end if
716 
717 !    ###########################################################
718 !    ### 13. Write the history into the _HIST file
719 !    ###
720 
721 #if defined HAVE_NETCDF
722      if (need_writeHIST.and.me==master) then
723        ifirst=merge(0,1,(itime>1.or.icycle>1))
724        call write_md_hist(hist,filename,ifirst,itime_hist,ab_mover%natom,scfcv_args%dtset%nctime,&
725 &       ab_mover%ntypat,ab_mover%typat,amu_curr,ab_mover%znucl,ab_mover%dtion,scfcv_args%dtset%mdtemp)
726      end if
727 #endif
728 
729 !    ###########################################################
730 !    ### 14. Output after SCFCV
731      if(need_verbose.and.need_scfcv_cycle)then
732        write(message,'(a,3a,a,72a)')&
733 &       ch10,('-',kk=1,3),'OUTPUT',('-',kk=1,71)
734        call wrtout(ab_out,message,'COLL')
735        call wrtout(std_out,message,'COLL')
736      end if
737      if (useprtxfase) then
738        call prtxfase(ab_mover,hist,itime_hist,ab_out,mover_AFTER)
739        call prtxfase(ab_mover,hist,itime_hist,std_out,mover_AFTER)
740      end if
741 
742 !    ###########################################################
743 !    ### 15. => Test Convergence of forces and stresses
744 
745      if (itime==ntime.and.icycle==ncycle)then
746        iexit=1
747        stat4xml="Failed"
748      else
749        stat4xml="Succeded"
750      end if
751 
752 !    Only if convergence is needed
753      if(specs%isFconv)then
754        if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5).or.mod(itime,2)==1)then
755          if (scfcv_args%dtset%tolmxf/=0)then
756            call fconv(hist%fcart(:,:,hist%ihist),&
757 &           scfcv_args%dtset%iatfix, &
758 &           iexit,itime,&
759 &           ab_mover%natom,&
760 &           ntime,&
761 &           ab_mover%optcell,&
762 &           scfcv_args%dtset%strfact,&
763 &           scfcv_args%dtset%strtarget,&
764 &           hist%strten(:,hist%ihist),&
765 &           scfcv_args%dtset%tolmxf)
766          else
767            call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde)
768          end if
769        end if
770      end if
771 
772      if (itime==ntime.and.icycle==ncycle) iexit=1
773 
774 !    ###########################################################
775 !    ### 16. => Precondition forces, stress and energy
776 !    ### 17. => Call to each predictor
777 
778      call precpred_1geo(ab_mover,ab_xfh,amu_curr,deloc,&
779 &     scfcv_args%dtset%chkdilatmx,&
780 &     scfcv_args%mpi_enreg%comm_cell,&
781 &     scfcv_args%dtset%dilatmx,dtfil%filnam_ds(4),&
782 &     hist,scfcv_args%dtset%hmctt,&
783 &     icycle,iexit,itime,mttk_vars,&
784 &     scfcv_args%dtset%nctime,ncycle,nerr_dilatmx,scfcv_args%dtset%npsp,ntime,&
785 &     scfcv_args%dtset%rprimd_orig,skipcycle,&
786 &     scfcv_args%dtset%usewvl)
787 
788 !    Write MOLDYN netcdf and POSABIN files (done every dtset%nctime time step)
789 !    This file is not created for multibinit run
790      if(need_scfcv_cycle .and. (ab_mover%ionmov/=23 .or. icycle==1))then
791        if (scfcv_args%dtset%nctime>0) then
792          jj=itime; if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1) jj=jj-hist_prev%mxhist
793          if (jj>0) then
794            option=3
795            ihist_prev = abihist_findIndex(hist,-1)
796            call wrt_moldyn_netcdf(ab_mover%amass,scfcv_args%dtset,jj,option,dtfil%fnameabo_moldyn,&
797 &           scfcv_args%mpi_enreg,scfcv_args%results_gs,&
798 &           hist%rprimd(:,:,ihist_prev),dtfil%unpos,hist%vel(:,:,hist%ihist),&
799 &           hist%xred(:,:,ihist_prev))
800          end if
801          if (iexit==1) hist%ihist=ihist_prev
802        end if
803      end if
804      if(iexit/=0) exit
805 
806 !    ###########################################################
807 !    ### 18. Use the history  to extract the new values
808 !    ###     acell, rprimd and xred
809 
810      call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
811 
812      if(ab_mover%optcell/=0)then
813 
814        call matr3inv(rprimd,gprimd)
815 
816 !      If metric has changed since the initialization, update the Ylm's
817        if (scfcv_args%psps%useylm==1)then
818          option=0;
819          if (scfcv_args%dtset%iscf>0) option=1
820          call initylmg(gprimd,&
821 &         scfcv_args%kg,&
822 &         scfcv_args%dtset%kptns,&
823 &         scfcv_args%dtset%mkmem,&
824 &         scfcv_args%mpi_enreg,&
825 &         scfcv_args%psps%mpsang,&
826 &         scfcv_args%dtset%mpw,&
827 &         scfcv_args%dtset%nband,&
828 &         scfcv_args%dtset%nkpt,&
829 &         scfcv_args%npwarr,&
830 &         scfcv_args%dtset%nsppol,&
831 &         option,rprimd,&
832 &         scfcv_args%ylm,&
833 &         scfcv_args%ylmgr)
834        end if
835 
836      end if
837 
838      vel(:,:)=hist%vel(:,:,hist%ihist)
839 
840 !    vel_cell(3,3)= velocities of cell parameters
841 !    Not yet used here but compute it for consistency
842      vel_cell(:,:)=zero
843      if (ab_mover%ionmov==13) then
844        if (itime_hist>2) then
845          ihist_prev2 = abihist_findIndex(hist,-2)
846          vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)- hist%rprimd(:,:,ihist_prev2))/(two*ab_mover%dtion)
847        else if (itime_hist>1) then
848          ihist_prev = abihist_findIndex(hist,-1)
849          vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)-hist%rprimd(:,:,ihist_prev))/(ab_mover%dtion)
850        end if
851      end if
852 
853 !    This is needed for some compilers such as
854 !    pathscale, g95, xlf that do not exit
855 !    from a loop if you change the upper limit
856 !    inside
857      if (icycle>=ncycle .and. scfcv_args%mpi_enreg%me == 0) then
858        if(need_verbose)write(std_out,*) 'EXIT:',icycle,ncycle
859        exit
860      end if
861 
862      write(message,*) 'ICYCLE',icycle,skipcycle
863      if(need_verbose)call wrtout(std_out,message,'COLL')
864      write(message,*) 'NCYCLE',ncycle
865      if(need_verbose)call wrtout(std_out,message,'COLL')
866      if (skipcycle) exit
867 
868 !DEBUG
869      write(std_out,*)' m_mover : will call precpred_1geo'
870 !    call flush(std_out)
871 !ENDDEBUG
872 
873 
874 !    ###########################################################
875 !    ### 19. End loop icycle
876 
877    end do ! do icycle=1,ncycle
878 
879    if(iexit/=0)exit
880 
881 !  ###########################################################
882 !  ### 20. End loop itime
883 
884  end do ! do itime=1,ntime
885 
886  ! Call fconv here if we exited due to wall time limit.
887  if (timelimit_exit==1 .and. specs%isFconv) then
888    iexit = timelimit_exit
889    ntime = itime-1
890    ihist_prev = abihist_findIndex(hist,-1)
891    if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5)) then
892      if (scfcv_args%dtset%tolmxf/=0)then
893        call fconv(hist%fcart(:,:,ihist_prev),&
894 &       scfcv_args%dtset%iatfix, &
895 &       iexit, itime,&
896 &       ab_mover%natom,&
897 &       ntime,&
898 &       ab_mover%optcell,&
899 &       scfcv_args%dtset%strfact,&
900 &       scfcv_args%dtset%strtarget,&
901 &       hist%strten(:,ihist_prev),&
902 &       scfcv_args%dtset%tolmxf)
903      else
904        call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde)
905      end if
906    end if
907  end if
908 
909  ! Avoid pending requests if itime == ntime.
910  call xmpi_wait(quitsum_request,ierr)
911 
912 !###########################################################
913 !### 21. Set the final values of xred with the last
914 !###     computed values (not the last predicted)
915 
916  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
917  vel(:,:)=hist%vel(:,:,hist%ihist)
918 
919  if (DEBUG .and. ab_mover%ionmov==1)then
920    write (std_out,*) 'vel'
921    do kk=1,ab_mover%natom
922      write (std_out,*) hist%vel(:,kk,hist%ihist)
923    end do
924  end if
925 
926 !###########################################################
927 !### 22. XML Output at the end
928 
929 !XML output of the status
930  if (scfcv_args%mpi_enreg%me == 0 .and. scfcv_args%dtset%prtxml == 1) then
931    write(ab_xml_out, "(3a)") '    <geometryMinimisation type="',trim(specs%type4xml),'">'
932    write(ab_xml_out, "(5a)") '      <status cvState="',trim(stat4xml) ,'" stop-criterion="',trim(specs%crit4xml),'" />'
933    write(ab_xml_out, "(3a)") '    </geometryMinimisation>'
934  end if
935 
936 !###########################################################
937 !### 23. Deallocate hist and ab_mover datatypes
938 
939 !This call is needed to free an internal matrix. However, this is not optimal ...
940 !One should instead have a datastructure associated with the preconditioner...
941  if (ab_mover%goprecon>0)then
942    call prec_simple(ab_mover,preconforstr,hist,1,1,1)
943  end if
944 
945  if (ab_mover%ionmov==13)then
946    call mttk_fin(mttk_vars)
947  end if
948 
949  if (ab_mover%ionmov==10 .or. ab_mover%ionmov==11)then
950    call delocint_fin(deloc)
951  end if
952 
953  ABI_DEALLOCATE(xred_prev)
954 
955  call abihist_free(hist)
956  call abihist_free(hist_prev)
957 
958  call abimover_destroy(ab_mover)
959  call abiforstr_fin(preconforstr)
960 
961  call status(0,dtfil%filstat,iexit,level,'exit          ')
962 
963 contains

ABINIT/prtnatom [ Functions ]

[ Top ] [ Functions ]

NAME

 prtnatom

FUNCTION

 Print information for N atoms

INPUTS

 prtallatoms = Logical for PRTint ALL ATOMS
 atlist      = ATom LIST
 index       = index for each atom
 natom       = Number of ATOMs

OUTPUT

  tag = The string to put for aech atom

PARENTS

      prtxfase

CHILDREN

      gettag,wrtout

SOURCE

1594 subroutine prtnatom(atlist,iout,message,natom,prtallatoms,thearray)
1595 
1596 
1597 !This section has been created automatically by the script Abilint (TD).
1598 !Do not modify the following lines by hand.
1599 #undef ABI_FUNC
1600 #define ABI_FUNC 'prtnatom'
1601 !End of the abilint section
1602 
1603 implicit none
1604 
1605 !Arguments ------------------------------------
1606 !scalars
1607   logical,intent(in) :: prtallatoms
1608   integer,intent(in) :: natom
1609   logical,intent(in) :: atlist(natom)
1610   integer,intent(in) :: iout
1611   character(len=*),intent(inout) :: message
1612 !arrays
1613   real(dp) :: thearray(3,natom)
1614 
1615 !Local variables-------------------------------
1616 !scalars
1617   integer :: kk
1618   character(len=7)   :: tag ! Maximal ' (9999)'
1619   character(len=18)   :: fmt
1620 
1621 ! *********************************************************************
1622 
1623  fmt='(a,a,1p,3e22.14,a)'
1624 
1625  do kk=1,natom
1626    if (atlist(kk)) then
1627      call gettag(atlist,kk,natom,prtallatoms,tag)
1628      write(message,fmt)TRIM(message),ch10,thearray(:,kk),tag
1629    end if
1630  end do
1631  !MGNAG Runtime Error: wrtout_cpp.f90, line 896: Buffer overflow on output
1632  call wrtout(iout,message,'COLL')
1633 
1634  end subroutine prtnatom

ABINIT/prtxfase [ Functions ]

[ Top ] [ Functions ]

NAME

 prtxfase

FUNCTION

 Print the values of xcart (X), forces (F) acell (A), Stresses (S), and energy (E)
 All values come from the history hist
 Also compute and print max and rms forces.
 Also compute absolute and relative differences with previous calculation

INPUTS

 ab_mover<type abimover>=Subset of dtset only related with
          |                 movement of ions and acell, contains:
          | dtion:  Time step
          ! natom:  Number of atoms
          | vis:    viscosity
          | iatfix: Index of atoms and directions fixed
          | amass:  Mass of ions
 hist<type abihist>=Historical record of positions, forces,
                               stresses, cell and energies,
 itime= time step
 iout=unit number for printing

OUTPUT

  (only writing)

PARENTS

      mover

CHILDREN

      gettag,wrtout

SOURCE

1253 subroutine prtxfase(ab_mover,hist,itime,iout,pos)
1254 
1255 
1256 !This section has been created automatically by the script Abilint (TD).
1257 !Do not modify the following lines by hand.
1258 #undef ABI_FUNC
1259 #define ABI_FUNC 'prtxfase'
1260 !End of the abilint section
1261 
1262 implicit none
1263 
1264 !Arguments ------------------------------------
1265 !scalars
1266  type(abimover),intent(in) :: ab_mover
1267  type(abihist),intent(in),target :: hist
1268  integer,intent(in) :: itime,iout
1269  integer,intent(in) :: pos
1270 !arrays
1271 
1272 !Local variables-------------------------------
1273 !scalars
1274  integer :: jj,kk,unfixd,iprt
1275  real(dp) :: val_max,val_rms,ucvol ! Values maximal and RMS, Volume of Unitary cell
1276  real(dp) :: dEabs,dErel ! Diff of energy absolute and relative
1277  real(dp) :: ekin
1278  real(dp) :: angle(3),rmet(3,3)
1279 !character(len=80*(max(ab_mover%natom,3)+1)) :: msg
1280 !MGNAG: This is not very safe. One should use line-based output istead of appending chars
1281 ! and then outputting everything! For the time being I use this temporary hack to solve the problem with NAG
1282  character(len=max(80*(max(ab_mover%natom,3)+1),50000)) :: msg
1283  character(len=18)   :: fmt1
1284  logical :: prtallatoms
1285 !arrays
1286  logical :: atlist(ab_mover%natom)
1287  real(dp),allocatable :: fred(:,:),xcart(:,:)
1288  real(dp),pointer :: acell(:),fcart(:,:),rprimd(:,:),strten(:),vel(:,:),xred(:,:)
1289 
1290 ! ***********************************************************
1291 
1292  fmt1='(a,a,1p,3e22.14)'
1293 
1294 !##########################################################
1295 !### 1. Organize list of atoms to print
1296 
1297  prtallatoms=.TRUE.
1298  do kk=1,ab_mover%natom
1299    if (ab_mover%prtatlist(kk)/=kk) prtallatoms=.FALSE.
1300  end do
1301 
1302  atlist(:)=.FALSE.
1303  do iprt=1,ab_mover%natom
1304    if (ab_mover%prtatlist(iprt)>0.and.ab_mover%prtatlist(iprt)<=ab_mover%natom) atlist(ab_mover%prtatlist(iprt))=.TRUE.
1305  end do
1306 
1307  acell  => hist%acell(:,hist%ihist)
1308  rprimd => hist%rprimd(:,:,hist%ihist)
1309  xred   => hist%xred(:,:,hist%ihist)
1310  fcart  => hist%fcart(:,:,hist%ihist)
1311  strten => hist%strten(:,hist%ihist)
1312  vel    => hist%vel(:,:,hist%ihist)
1313 
1314 !###########################################################
1315 !### 1. Positions
1316 
1317  ABI_ALLOCATE(xcart,(3,ab_mover%natom))
1318  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
1319 
1320  write(msg, '(a,a)' )ch10,' Cartesian coordinates (xcart) [bohr]'
1321  call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xcart)
1322 
1323  write(msg, '(a)' )' Reduced coordinates (xred)'
1324  call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xred)
1325 
1326  ABI_DEALLOCATE(xcart)
1327 
1328 !###########################################################
1329 !### 2. Forces
1330 
1331  if(pos==mover_AFTER)then
1332 
1333    ABI_ALLOCATE(fred,(3,ab_mover%natom))
1334    call fcart2fred(fcart,fred,rprimd,ab_mover%natom)
1335 
1336 !  Compute max |f| and rms f,
1337 !  EXCLUDING the components determined by iatfix
1338    val_max=0.0_dp
1339    val_rms=0.0_dp
1340    unfixd=0
1341    do kk=1,ab_mover%natom
1342      do jj=1,3
1343        if (ab_mover%iatfix(jj,kk) /= 1) then
1344          unfixd=unfixd+1
1345          val_rms=val_rms+fcart(jj,kk)**2
1346          val_max=max(val_max,abs(fcart(jj,kk)**2))
1347        end if
1348      end do
1349    end do
1350    if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd))
1351 
1352    write(msg, '(a,1p,2e12.5,a)' ) ' Cartesian forces (fcart) [Ha/bohr]; max,rms=',sqrt(val_max),val_rms,' (free atoms)'
1353    call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,fcart)
1354 
1355    write(msg, '(a)' )' Reduced forces (fred)'
1356    call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,fred)
1357    ABI_DEALLOCATE(fred)
1358  end if
1359 
1360 !###########################################################
1361 !### 3. Velocities
1362 
1363 !Only if the velocities are being used
1364  if (hist%isVused)then
1365 !  Only if velocities are recorded in a history
1366    if (allocated(hist%vel))then
1367 !    Compute max |v| and rms v,
1368 !    EXCLUDING the components determined by iatfix
1369      val_max=0.0_dp
1370      val_rms=0.0_dp
1371      unfixd=0
1372      do kk=1,ab_mover%natom
1373        do jj=1,3
1374          if (ab_mover%iatfix(jj,kk) /= 1) then
1375            unfixd=unfixd+1
1376            val_rms=val_rms+vel(jj,kk)**2
1377            val_max=max(val_max,abs(vel(jj,kk)**2))
1378          end if
1379        end do
1380      end do
1381      if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd))
1382 
1383      write(msg, '(a,1p,2e12.5,a)' ) &
1384 &     ' Cartesian velocities (vel) [bohr*Ha/hbar]; max,rms=',sqrt(val_max),val_rms,' (free atoms)'
1385      call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,vel)
1386 
1387 !    Compute the ionic kinetic energy (no cell shape kinetic energy yet)
1388      ekin=0.0_dp
1389      do kk=1,ab_mover%natom
1390        do jj=1,3
1391 !        Warning : the fixing of atoms is implemented in reduced
1392 !        coordinates, so that this expression is wrong
1393          if (ab_mover%iatfix(jj,kk) == 0) then
1394            ekin=ekin+0.5_dp*ab_mover%amass(kk)*vel(jj,kk)**2
1395          end if
1396        end do
1397      end do
1398      write(msg, '(a,1p,e22.14,a)' )' Kinetic energy of ions (ekin) [Ha]=',ekin
1399      call wrtout(iout,msg,'COLL')
1400    end if
1401  end if
1402 
1403 !###########################################################
1404 !### 3. ACELL
1405 
1406 !Only if the acell is being used
1407  if (hist%isARused)then
1408 !  Only if acell is recorded in a history
1409    if (allocated(hist%acell))then
1410      write(msg, '(a)' )' Scale of Primitive Cell (acell) [bohr]'
1411      write(msg,fmt1)TRIM(msg),ch10,acell(:)
1412      call wrtout(iout,msg,'COLL')
1413    end if
1414  end if
1415 
1416 !###########################################################
1417 !### 4. RPRIMD
1418 
1419 !Only if the acell is being used
1420  if (hist%isARused)then
1421 !  Only if rprimd is recorded in a history
1422    if (allocated(hist%rprimd))then
1423      write(msg, '(a)' )' Real space primitive translations (rprimd) [bohr]'
1424      do kk=1,3
1425        write(msg,fmt1)TRIM(msg),ch10,rprimd(:,kk)
1426      end do
1427      call wrtout(iout,msg,'COLL')
1428    end if
1429  end if
1430 
1431 !###########################################################
1432 !### 5. Unitary cell volume
1433 
1434  if (ab_mover%optcell/=0)then
1435 
1436    ucvol=&
1437 &   rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
1438 &   rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
1439 &   rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
1440 
1441    write(msg, '(a,1p,e22.14)' )' Unitary Cell Volume (ucvol) [Bohr^3]=',ucvol
1442    call wrtout(iout,msg,'COLL')
1443 
1444 !  ###########################################################
1445 !  ### 5. Angles and lengths
1446 
1447 !  Compute real space metric.
1448    rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
1449 
1450    angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0
1451    angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0
1452    angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0
1453 
1454    write(msg, '(a)' )' Angles (23,13,12)= [degrees]'
1455    write(msg,fmt1)TRIM(msg),ch10,angle(:)
1456    call wrtout(iout,msg,'COLL')
1457 
1458    write(msg, '(a)' ) ' Lengths [Bohr]'
1459    write(msg,fmt1)TRIM(msg),ch10,sqrt(rmet(1,1)),sqrt(rmet(2,2)),sqrt(rmet(3,3))
1460    call wrtout(iout,msg,'COLL')
1461 
1462 !  ###########################################################
1463 !  ### 5. Stress Tensor
1464 
1465    if(pos==mover_AFTER)then
1466 !    Only if strten is recorded in a history
1467      if (allocated(hist%strten))then
1468 
1469        write(msg, '(a)' ) ' Stress tensor in cartesian coordinates (strten) [Ha/bohr^3]'
1470 
1471        write(msg,fmt1)TRIM(msg),ch10,strten(1),strten(6),strten(5)
1472        write(msg,fmt1)TRIM(msg),ch10,strten(6),strten(2),strten(4)
1473        write(msg,fmt1)TRIM(msg),ch10,strten(5),strten(4),strten(3)
1474        call wrtout(iout,msg,'COLL')
1475      end if
1476    end if
1477  end if
1478 
1479 !###########################################################
1480 !### 6. Energy
1481 
1482  if(pos==mover_AFTER)then
1483    write(msg, '(a,1p,e22.14)' )' Total energy (etotal) [Ha]=',hist%etot(hist%ihist)
1484 
1485    if (itime>1)then
1486      jj = abihist_findIndex(hist,-1)
1487      dEabs=hist%etot(hist%ihist)-hist%etot(jj)
1488      dErel=2*dEabs/(abs(hist%etot(hist%ihist))+abs(hist%etot(jj)))
1489      write(msg, '(a,a,a,a)' )TRIM(msg),ch10,ch10,' Difference of energy with previous step (new-old):'
1490      write(msg, '(a,a,10a,a,1p,e12.5,a,10a,a,1p,e12.5)')&
1491 &     TRIM(msg),ch10,&
1492 &     (' ',jj=1,10),' Absolute (Ha)=',dEabs,ch10,&
1493 &     (' ',jj=1,10),' Relative     =',dErel
1494    end if
1495    call wrtout(iout,msg,'COLL')
1496  end if
1497 
1498  contains

ABINIT/wrt_moldyn_netcdf [ Functions ]

[ Top ] [ Functions ]

NAME

 wrt_moldyn_netcdf

FUNCTION

 Write two files for later molecular dynamics analysis:
  - MOLDYN.nc (netcdf format) : evolution of key quantities with time (pressure, energy, ...)
  - POSABIN : values of coordinates and velocities for the next time step

COPYRIGHT

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

INPUTS

  amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
  dtset <type(dataset_type)>=all input variables for this dataset
  itime=time step index
  option=1: write MOLDYN.nc file (netcdf format)
         2: write POSABIN file
         3: write both
  moldyn_file=name of the MD netcdf file
  mpi_enreg=informations about MPI parallelization
  results_gs <type(results_gs_type)>=results (energy and its components,
   forces and its components, the stress tensor) of a ground-state computation
  rprimd(3,3)=real space primitive translations
  unpos=unit number for POSABIN file
  vel(3,natom)=velocities of atoms
  xred(3,natom)=reduced coordinates of atoms

OUTPUT

  -- only printing --

SIDE EFFECTS

PARENTS

      mover

CHILDREN

      metric,wrtout,xcart2xred,xred2xcart

SOURCE

1685 subroutine wrt_moldyn_netcdf(amass,dtset,itime,option,moldyn_file,mpi_enreg,&
1686 &                            results_gs,rprimd,unpos,vel,xred)
1687 
1688  use defs_basis
1689  use defs_abitypes
1690  use m_results_gs
1691  use m_abicore
1692  use m_errors
1693 #if defined HAVE_NETCDF
1694  use netcdf
1695 #endif
1696 
1697  use m_io_tools,   only : open_file
1698  use m_geometry,   only : xcart2xred, xred2xcart, metric
1699 
1700 !This section has been created automatically by the script Abilint (TD).
1701 !Do not modify the following lines by hand.
1702 #undef ABI_FUNC
1703 #define ABI_FUNC 'wrt_moldyn_netcdf'
1704 !End of the abilint section
1705 
1706  implicit none
1707 
1708 !Arguments ------------------------------------
1709 !scalars
1710  integer,intent(in) :: itime,option,unpos
1711  character(fnlen),intent(in) :: moldyn_file
1712  type(dataset_type),intent(in) :: dtset
1713  type(MPI_type),intent(in) :: mpi_enreg
1714  type(results_gs_type),intent(in) :: results_gs
1715 !arrays
1716  real(dp),intent(in) :: amass(dtset%natom),rprimd(3,3)
1717  real(dp),intent(in),target :: vel(3,dtset%natom)
1718  real(dp),intent(in) :: xred(3,dtset%natom)
1719 
1720 !Local variables-------------------------------
1721 !scalars
1722  integer,save :: ipos=0
1723  integer :: iatom,ii
1724  character(len=500) :: message
1725 #if defined HAVE_NETCDF
1726  integer :: AtomNumDimid,AtomNumId,CelId,CellVolumeId,DimCoordid,DimScalarid,DimVectorid
1727  integer :: EkinDimid,EkinId,EpotDimid,EpotId,EntropyDimid,EntropyId,MassDimid,MassId,NbAtomsid
1728  integer :: ncerr,ncid,PosId,StressDimid,StressId,TensorSymDimid
1729  integer :: TimeDimid,TimestepDimid,TimestepId
1730  logical :: atom_fix
1731  real(dp) :: ekin,ucvol
1732  character(len=fnlen) :: ficname
1733  character(len=16) :: chain
1734 #endif
1735 !arrays
1736 #if defined HAVE_NETCDF
1737  integer :: PrimVectId(3)
1738  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1739  real(dp),allocatable ::  xcart(:,:)
1740  real(dp),pointer :: vcart(:,:),vred(:,:),vtmp(:,:)
1741 #endif
1742 
1743 ! *************************************************************************
1744 
1745 !Only done by master processor, every nctime step
1746  if (mpi_enreg%me==0.and.dtset%nctime>0) then
1747 
1748 !  Netcdf file name
1749 #if defined HAVE_NETCDF
1750    ficname = trim(moldyn_file)//'.nc'
1751 #endif
1752 
1753 !  Xcart from Xred
1754 #if defined HAVE_NETCDF
1755    ABI_ALLOCATE(xcart,(3,dtset%natom))
1756    call xred2xcart(dtset%natom,rprimd,xcart,xred)
1757 #endif
1758 
1759 !  ==========================================================================
1760 !  First time step: write header of netcdf file
1761 !  ==========================================================================
1762    if (itime==1.and.(option==1.or.option==3)) then
1763 
1764      ipos=0
1765 
1766 #if defined HAVE_NETCDF
1767 !    Write message
1768      write(message,'(4a)')ch10,' Open file ',trim(ficname),' to store molecular dynamics information.'
1769      call wrtout(std_out,message,'COLL')
1770 
1771 !    Create netcdf file
1772      ncerr = nf90_create(ficname, NF90_CLOBBER , ncid)
1773      NCF_CHECK_MSG(ncerr,'nf90_create')
1774 
1775 !    Dimension time for netcdf (time dim is unlimited)
1776      ncerr = nf90_def_dim(ncid, "time", nf90_unlimited, TimeDimid)
1777      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1778 
1779 !    Symetric Tensor Dimension
1780      ncerr = nf90_def_dim(ncid, "DimTensor", size(results_gs%strten), TensorSymDimid)
1781      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1782 
1783 !    Coordinates Dimension
1784      ncerr = nf90_def_dim(ncid, "DimCoord", size(xcart,1), DimCoordid)
1785      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1786 
1787 !    Atoms Dimensions
1788      ncerr = nf90_def_dim(ncid, "NbAtoms", dtset%natom, NbAtomsid)
1789      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1790 
1791 !    Vector Dimension
1792      ncerr = nf90_def_dim(ncid, "DimVector", 3 , DimVectorid)
1793      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1794 
1795 !    Scalar Dimension
1796      ncerr = nf90_def_dim(ncid, "DimScalar", 1 , DimScalarid)
1797      NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1798 
1799 !    Time step and time unit
1800      ncerr = nf90_def_var(ncid, "Time_step", nf90_double , DimScalarid, TimestepDimid)
1801      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1802      ncerr = nf90_put_att(ncid, TimestepDimid, "units", "atomic time unit")
1803      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1804 
1805 !    Ionic masses
1806      ncerr = nf90_def_var(ncid, "Ionic_Mass", nf90_double , NbAtomsid, MassDimid)
1807      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1808      ncerr = nf90_put_att(ncid, MassDimid, "units", "atomic mass unit")
1809      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1810 
1811 !    Ionic atomic numbers
1812      ncerr = nf90_def_var(ncid, "Ionic_Atomic_Number", nf90_double , NbAtomsid, AtomNumDimid)
1813      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1814 
1815 !    E_pot
1816      ncerr = nf90_def_var(ncid, "E_pot", nf90_double , TimeDimid, EpotDimid)
1817      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1818      ncerr = nf90_put_att(ncid, EpotDimid, "units", "hartree")
1819      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1820 
1821 !    E_kin
1822      ncerr = nf90_def_var(ncid, "E_kin", nf90_double , TimeDimid, EkinDimid)
1823      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1824      ncerr = nf90_put_att(ncid, EkinDimid, "units", "hartree")
1825      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1826 
1827 !    Entropy
1828      ncerr = nf90_def_var(ncid, "Entropy", nf90_double , TimeDimid, EntropyDimid)
1829      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1830      ncerr = nf90_put_att(ncid, EntropyDimid, "units", "")
1831      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1832 
1833 !    Stress tensor
1834      ncerr = nf90_def_var(ncid, "Stress", nf90_double , (/TensorSymDimid,TimeDimid/), StressDimid)
1835      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1836      ncerr = nf90_put_att(ncid, StressDimid, "units", "hartree/bohr^3")
1837      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1838 
1839 !    Positions
1840      ncerr = nf90_def_var(ncid, "Position", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), PosId)
1841      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1842      ncerr = nf90_put_att(ncid, PosId, "units", "bohr")
1843      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1844 
1845 !    Celerities
1846      ncerr = nf90_def_var(ncid, "Celerity", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), CelId)
1847      NCF_CHECK_MSG(ncerr,'nf90_def_var')
1848      ncerr = nf90_put_att(ncid, CelId, "units", "bohr/(atomic time unit)")
1849      NCF_CHECK_MSG(ncerr,'nf90_put_att')
1850 
1851 !    In case of volume cell constant
1852      if (dtset%optcell==0) then
1853 !      Primitive vectors
1854        do ii = 1,3
1855          write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii
1856          ncerr = nf90_def_var(ncid, trim(chain), nf90_double , DimVectorid, PrimVectId(ii))
1857          NCF_CHECK_MSG(ncerr,'nf90_def_var')
1858        end do
1859 !      Cell Volume
1860        ncerr = nf90_def_var(ncid, "Cell_Volume", nf90_double , DimScalarid, CellVolumeId)
1861        NCF_CHECK_MSG(ncerr,'nf90_def_var')
1862        ncerr = nf90_put_att(ncid, CellVolumeId, "units", "bohr^3")
1863        NCF_CHECK_MSG(ncerr,'nf90_put_att')
1864      end if
1865 
1866 !    Leave define mode and close file
1867      ncerr = nf90_enddef(ncid)
1868      NCF_CHECK_MSG(ncerr,'nf90_enddef')
1869      ncerr = nf90_close(ncid)
1870      NCF_CHECK_MSG(ncerr,'nf90_close')
1871 #endif
1872    end if
1873 
1874 !  ==========================================================================
1875 !  Write data to netcdf file (every nctime time step)
1876 !  ==========================================================================
1877    if (mod(itime, dtset%nctime)==0.and.(option==1.or.option==3)) then
1878 
1879      ipos=ipos+1
1880 
1881 #if defined HAVE_NETCDF
1882 !    Write message
1883      write(message,'(3a)')ch10,' Store molecular dynamics information in file ',trim(ficname)
1884      call wrtout(std_out,message,'COLL')
1885 
1886 !    Open netcdf file
1887      ncerr = nf90_open(ficname, nf90_write, ncid)
1888      NCF_CHECK_MSG(ncerr,'nf90_open')
1889 
1890 !    Time step
1891      ncerr = nf90_inq_varid(ncid, "Time_step", TimestepId)
1892      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1893      ncerr = nf90_put_var(ncid, TimestepId, dtset%dtion)
1894      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1895 
1896 !    Ionic masses
1897      ncerr = nf90_inq_varid(ncid, "Ionic_Mass", MassId)
1898      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1899      ncerr = nf90_put_var(ncid, MassId, amass, start = (/ 1 /), count=(/dtset%natom/))
1900      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1901 
1902 !    Ionic atomic numbers
1903      ncerr = nf90_inq_varid(ncid, "Ionic_Atomic_Number", AtomNumId)
1904      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1905      ncerr = nf90_put_var(ncid, AtomNumId, dtset%znucl(dtset%typat(:)),start=(/1/),count=(/dtset%natom/))
1906      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1907 
1908 !    Epot
1909      ncerr = nf90_inq_varid(ncid, "E_pot", EpotId)
1910      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1911      ncerr = nf90_put_var(ncid, EpotId, (/results_gs%etotal/), start=(/ipos/),count=(/1/))
1912      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1913 
1914 !    Ekin
1915      ekin=zero;atom_fix=(maxval(dtset%iatfix)>0)
1916      if (dtset%ionmov==1.or.(.not.atom_fix)) then
1917        vcart => vel
1918      else
1919        ABI_ALLOCATE(vcart,(3,dtset%natom))
1920        ABI_ALLOCATE(vred,(3,dtset%natom))
1921        vtmp => vel
1922        call xcart2xred(dtset%natom,rprimd,vtmp,vred)
1923        do iatom=1,dtset%natom
1924          do ii=1,3
1925            if (dtset%iatfix(ii,iatom)==1) vred(ii,iatom)=zero
1926          end do
1927        end do
1928        call xred2xcart(dtset%natom,rprimd,vcart,vred)
1929        ABI_DEALLOCATE(vred)
1930      end if
1931      do iatom=1,dtset%natom
1932        do ii=1,3
1933          ekin=ekin+half*amass(iatom)*vcart(ii,iatom)**2
1934        end do
1935      end do
1936      if (dtset%ionmov/=1.and.atom_fix)  then
1937        ABI_DEALLOCATE(vcart)
1938      end if
1939      ncerr = nf90_inq_varid(ncid, "E_kin", EkinId)
1940      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1941      ncerr = nf90_put_var(ncid, EkinId, (/ekin/), start = (/ipos/),count=(/1/))
1942      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1943 
1944 !    EntropyDimid
1945      ncerr = nf90_inq_varid(ncid, "Entropy", EntropyId)
1946      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1947      ncerr = nf90_put_var(ncid, EntropyId, (/results_gs%energies%entropy/),start = (/ipos/),count=(/1/))
1948      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1949 
1950 !    Stress tensor
1951      ncerr = nf90_inq_varid(ncid, "Stress", StressId)
1952      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1953      ncerr = nf90_put_var(ncid, StressId, results_gs%strten, &
1954 &     start=(/1,ipos/),count=(/size(results_gs%strten)/))
1955      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1956 
1957 !    Positions
1958      ncerr = nf90_inq_varid(ncid, "Position", PosId)
1959      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1960      ncerr = nf90_put_var(ncid, PosId, xcart, start=(/1,1,ipos/), &
1961 &     count=(/size(xcart,1),dtset%natom,1/))
1962      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1963 
1964 !    Celerities
1965      ncerr = nf90_inq_varid(ncid, "Celerity", CelId)
1966      NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1967      ncerr = nf90_put_var(ncid, CelId, vel, start=(/1,1,ipos/), &
1968 &     count=(/size(vel,1),dtset%natom,1/)  )
1969      NCF_CHECK_MSG(ncerr,'nf90_put_var')
1970 
1971 !    In case of volume cell constant
1972      if (dtset%optcell==0.and.ipos==1) then
1973 !      Primitive vectors
1974        do ii = 1,3
1975          write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii
1976          ncerr = nf90_inq_varid(ncid, trim(chain), PrimVectId(ii) )
1977          NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1978          ncerr = nf90_put_var(ncid, PrimVectId(ii), rprimd(:,ii))
1979          NCF_CHECK_MSG(ncerr,'nf90_put_var')
1980        end do
1981 !      Cell Volume
1982        call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1983        ncerr = nf90_inq_varid(ncid, "Cell_Volume" , CellVolumeId)
1984        NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1985        ncerr = nf90_put_var(ncid, CellVolumeId, ucvol)
1986        NCF_CHECK_MSG(ncerr,'nf90_put_var')
1987      end if
1988 
1989 !    Close file
1990      ncerr = nf90_close(ncid)
1991      NCF_CHECK_MSG(ncerr,'nf90_close')
1992 #endif
1993    end if
1994 
1995 !  ==========================================================================
1996 !  Write data to POSABIN file (every nctime time step if option=3)
1997 !  ==========================================================================
1998    if ((mod(itime, dtset%nctime)==0.and.option==3).or.(option==2)) then
1999 
2000 !    Open file for writing
2001      if (open_file('POSABIN',message,unit=unpos,status='replace',form='formatted') /= 0 ) then
2002        MSG_ERROR(message)
2003      end if
2004 
2005 !    Write Positions
2006      if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'xred  ',(xred(ii,1),ii=1,3)
2007      if (dtset%natom>1) then
2008        do iatom=2,dtset%natom
2009          write(unpos,'(7x,3d18.5)') (xred(ii,iatom),ii=1,3)
2010        end do
2011      end if
2012 
2013 !    Write Velocities
2014      if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'vel  ',(vel(ii,1),ii=1,3)
2015      if (dtset%natom>1) then
2016        do iatom=2,dtset%natom
2017          write(unpos,'(7x,3d18.5)') (vel(ii,iatom),ii=1,3)
2018        end do
2019      end if
2020 
2021 !    Close file
2022      close(unpos)
2023    end if
2024 
2025 #if defined HAVE_NETCDF
2026    ABI_DEALLOCATE(xcart)
2027 #endif
2028 
2029 !  ==========================================================================
2030 !  End if master proc
2031  end if
2032 
2033 !Fake lines
2034 #if !defined HAVE_NETCDF
2035  if (.false.) write(std_out,*) moldyn_file,results_gs%etotal,rprimd(1,1)
2036 #endif
2037 
2038 end subroutine wrt_moldyn_netcdf