TABLE OF CONTENTS


ABINIT/anaddb [ Programs ]

[ Top ] [ Programs ]

NAME

 anaddb

FUNCTION

 Main routine for analysis of the interatomic force constants and associated properties.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (XG,DCA,JCC,CL,XW,GA)
 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

  (main routine)

OUTPUT

  (main routine)

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,anaddb_dtset_free,anaddb_init
      asrq0_apply,asrq0_free,crystal_free,ddb_diel,ddb_elast,ddb_free
      ddb_from_file,ddb_hdr_free,ddb_hdr_open_read,ddb_internalstr
      ddb_interpolate,ddb_piezo,dfpt_phfrq,dfpt_prtph,dfpt_symph
      elast_ncwrite,electrooptic,elphon,flush_unit,gruns_anaddb,gtblk9,gtdyn9
      harmonic_thermo,herald,ifc_free,ifc_init,ifc_outphbtrap,ifc_print
      ifc_speedofsound,ifc_write,instrng,int2char4,inupper,invars9,isfile
      mkphbs,mkphdos,nctk_defwrite_nonana_terms,outvars_anaddb,phdos_free
      phdos_ncwrite,phdos_print,phdos_print_debye,phdos_print_msqd
      phdos_print_thermo,ramansus,relaxpol,thermal_supercell_free
      thermal_supercell_make,thermal_supercell_print,thmeig,timein,wrtout
      xmpi_bcast,xmpi_end,xmpi_init,xmpi_sum

SOURCE

 40 #if defined HAVE_CONFIG_H
 41 #include "config.h"
 42 #endif
 43 
 44 #include "abi_common.h"
 45 
 46 program anaddb
 47 
 48  use defs_basis
 49  use m_build_info
 50  use m_xmpi
 51  use m_xomp
 52  use m_abicore
 53  use m_errors
 54  !use m_argparse
 55  use m_ifc
 56  use m_ddb
 57  use m_ddb_hdr
 58  use m_phonons
 59  use m_gruneisen
 60  use m_supercell
 61  use iso_c_binding
 62  use m_nctk
 63 #ifdef HAVE_NETCDF
 64  use netcdf
 65 #endif
 66 
 67  use m_io_tools,       only : open_file, flush_unit
 68  use m_fstrings,       only : int2char4, itoa, sjoin, strcat, inupper
 69  use m_specialmsg,     only : specialmsg_getcount, herald
 70  use m_time,           only : asctime, timein
 71  use m_parser,         only : instrng
 72  use m_dtfil,          only : isfile
 73  use m_anaddb_dataset, only : anaddb_init, anaddb_dataset_type, anaddb_dtset_free, outvars_anaddb, invars9
 74  use m_ddb_interpolate, only : ddb_interpolate
 75  use m_crystal,        only : crystal_t, crystal_free
 76  use m_crystal_io,     only : crystal_ncwrite
 77  use m_dynmat,         only : gtdyn9, dfpt_phfrq, dfpt_prtph
 78  use m_elphon,         only : elphon
 79  use m_harmonic_thermo,only : harmonic_thermo
 80  use m_thmeig,         only : thmeig
 81  use m_raman,          only : ramansus, electrooptic
 82  use m_ddb_diel,       only : ddb_diel
 83  use m_relaxpol,       only : relaxpol
 84  use m_ddb_elast,      only : ddb_elast
 85  use m_ddb_piezo,      only : ddb_piezo
 86  use m_ddb_internalstr, only : ddb_internalstr
 87 
 88 !This section has been created automatically by the script Abilint (TD).
 89 !Do not modify the following lines by hand.
 90 #undef ABI_FUNC
 91 #define ABI_FUNC 'anaddb'
 92 !End of the abilint section
 93 
 94  implicit none
 95 
 96 !Local variables-------------------------------
 97  integer :: msym !  msym =maximum number of symmetry elements in space group
 98 !Define input and output unit numbers (some are defined in defs_basis -all should be there ...):
 99  integer,parameter :: ddbun=2,master=0 ! FIXME: these should not be reserved unit numbers!
100  integer,parameter :: rftyp4=4
101  integer :: comm,iatom,iblok,iblok_stress,iblok_tmp,idir,ii,index
102  integer :: ierr,iphl2,lenstr,mtyp,mpert,msize,natom
103  integer :: nsym,ntypat,option,usepaw,nproc,my_rank,ana_ncid,prt_internalstr
104  logical :: iam_master
105  integer :: rfelfd(4),rfphon(4),rfstrs(4),ngqpt_coarse(3)
106  integer :: count_wminmax(2)
107  integer,allocatable :: d2flg(:)
108  real(dp) :: etotal,tcpu,tcpui,twall,twalli
109  real(dp) :: dielt(3,3)
110  real(dp) :: compl(6,6),compl_clamped(6,6),compl_stress(6,6)
111  real(dp) :: dielt_rlx(3,3),elast(6,6),elast_clamped(6,6),elast_stress(6,6)
112  real(dp) :: epsinf(3,3),red_ptot(3),pel(3)
113  real(dp) :: piezo(6,3),qphnrm(3),qphon(3,3),strten(6),tsec(2)
114  real(dp) :: wminmax(2)
115  real(dp),allocatable :: d2cart(:,:),dchide(:,:,:)
116  real(dp),allocatable :: dchidt(:,:,:,:),displ(:),eigval(:,:)
117  real(dp),allocatable :: eigvec(:,:,:,:,:),fact_oscstr(:,:,:),instrain(:,:)
118  real(dp),allocatable :: fred(:,:),lst(:),phfrq(:)
119  real(dp),allocatable :: rsus(:,:,:)
120  real(dp),allocatable :: zeff(:,:,:)
121  character(len=10) :: procstr
122  character(len=24) :: codename
123  character(len=24) :: start_datetime
124  character(len=strlen) :: string
125  character(len=fnlen) :: filnam(7),elph_base_name,tmpfilename
126  character(len=500) :: message
127  !type(args_t) :: args
128  type(anaddb_dataset_type) :: inp
129  type(phonon_dos_type) :: Phdos
130  type(ifc_type) :: Ifc,Ifc_coarse
131  type(ddb_type) :: ddb
132  type(ddb_hdr_type) :: ddb_hdr
133  type(asrq0_t) :: asrq0
134  type(crystal_t) :: Crystal
135  type(supercell_type), allocatable :: thm_scells(:)
136 #ifdef HAVE_NETCDF
137  integer :: phdos_ncid, ncerr
138 #endif
139 
140 !******************************************************************
141 
142  ! Change communicator for I/O (mandatory!)
143  call abi_io_redirect(new_io_comm=xmpi_world)
144 
145  ! Initialize MPI
146  call xmpi_init()
147 
148  ! MPI variables
149  comm = xmpi_world; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
150  iam_master = (my_rank == master)
151 
152  ! TODO
153  ! Parse command line arguments.
154  !args = args_parser(); if (args%exit /= 0) goto 100
155 
156 !Initialize memory profiling if it is activated !if a full abimem.mocc report is desired,
157 !set the argument of abimem_init to "2" instead of "0"
158 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
159 #ifdef HAVE_MEM_PROFILING
160  call abimem_init(0)
161 #endif
162 
163  ! Initialisation of the timing
164  call timein(tcpui,twalli)
165 
166  if (iam_master) then
167    codename='ANADDB'//repeat(' ',18)
168    call herald(codename,abinit_version,std_out)
169  end if
170 
171  start_datetime = asctime()
172 
173  ! Initialise the code: write heading, and read names of files.
174  if (iam_master) call anaddb_init(filnam)
175  call xmpi_bcast (filnam, master, comm, ierr)
176 
177  ! make log file for non-master procs
178  if (.not. iam_master) then
179    call int2char4(my_rank, procstr)
180    ABI_CHECK((procstr(1:1)/='#'),'Bug: string length too short!')
181    tmpfilename = trim(filnam(2)) // "_LOG_P" // trim(procstr)
182    if (open_file(tmpfilename, message, unit=std_out, form="formatted", action="write") /= 0) then
183      MSG_ERROR(message)
184    end if
185  end if
186 
187 !******************************************************************
188 
189  ! Must read natom from the DDB before being able to allocate some arrays needed for invars9
190 
191  call ddb_hdr_open_read(ddb_hdr,filnam(3),ddbun,DDB_VERSION,comm=comm, &
192 & dimonly=1)
193 
194  natom = ddb_hdr%natom
195  ntypat = ddb_hdr%ntypat
196  mtyp = ddb_hdr%mblktyp
197  usepaw = ddb_hdr%usepaw
198 
199  call ddb_hdr_free(ddb_hdr)
200 
201  mpert=natom+6
202  msize=3*mpert*3*mpert; if (mtyp==3) msize=msize*3*mpert
203 
204  ! Read the input file, and store the information in a long string of characters
205  ! strlen from defs_basis module
206  option=1
207  if (iam_master) then
208    call instrng (filnam(1),lenstr,option,strlen,string)
209 
210    ! To make case-insensitive, map characters to upper case.
211    call inupper(string(1:lenstr))
212  end if
213 
214  call xmpi_bcast(string,master, comm, ierr)
215  call xmpi_bcast(lenstr,master, comm, ierr)
216 
217  ! Read the inputs
218  call invars9(inp,lenstr,natom,string)
219 
220  !if (args%dry_run /= 0) then
221  !  call wrtout(std_out, "Dry run mode. Exiting after have read the input")
222  !  goto 100
223  !end if
224 
225  ! Open output files and ab_out (might change its name if needed)
226  ! MJV 1/2010 : now output file is open, but filnam(2) continues unmodified
227  ! so the other output files are overwritten instead of accumulating.
228  if (iam_master) then
229    tmpfilename = filnam(2)
230    call isfile(tmpfilename,'new')
231    if (open_file(tmpfilename,message,unit=ab_out,form='formatted',status='new') /= 0) then
232      MSG_ERROR(message)
233    end if
234    rewind (unit=ab_out)
235    call herald(codename,abinit_version,ab_out)
236 
237    ! Echo the inputs to console and main output file
238    call outvars_anaddb(inp,std_out)
239    call outvars_anaddb(inp,ab_out)
240  else
241    ab_out = dev_null
242  end if
243 
244 !******************************************************************
245 
246  ! Read the DDB information, also perform some checks, and symmetrize partially the DDB
247  write(message, '(a,a)' )' read the DDB information and perform some checks',ch10
248  call wrtout(std_out,message,'COLL')
249  call wrtout(ab_out,message,'COLL')
250 
251  ! DEBUG
252  !write(*,*) 'anaddb: natom=', natom
253  ! END DEBUG
254  call ddb_from_file(ddb,filnam(3),inp%brav,natom,inp%natifc,inp%atifc,Crystal,comm, prtvol=inp%prtvol)
255  nsym = Crystal%nsym
256  ! DEBUG
257  !do ii=1,Crystal%ntypat
258  !  write(*,*)'anaddb: amu=', crystal%amu(ii)
259  !end do
260  ! END DEBUG
261 
262  ! Acoustic Sum Rule
263  ! In case the interatomic forces are not calculated, the
264  ! ASR-correction (asrq0%d2asr) has to be determined here from the Dynamical matrix at Gamma.
265  asrq0 = ddb_get_asrq0(ddb, inp%asr, inp%rfmeth, crystal%xcart)
266 
267  ! TODO: This is to maintain the previous behaviour in which all the arrays were initialized to zero.
268  ! In the new version asrq0%d2asr is always computed if the Gamma block is present
269  ! and this causes changes in [v5][t28]
270  if (.not. (inp%ifcflag==0 .or. inp%instrflag/=0 .or. inp%elaflag/=0)) then
271    asrq0%d2asr = zero
272    if (asrq0%asr==3.or.asrq0%asr==4) then
273      asrq0%singular = zero; asrq0%uinvers = zero; asrq0%vtinvers = zero
274    end if
275  end if
276 
277  ! Open the netcdf file that will contain the anaddb results
278  ana_ncid = nctk_noid
279  if (iam_master) then
280 #ifdef HAVE_NETCDF
281    NCF_CHECK_MSG(nctk_open_create(ana_ncid, "anaddb.nc", xmpi_comm_self), "Creating anaddb.nc")
282    ncerr = nctk_def_dims(ana_ncid, [ &
283        nctkdim_t('number_of_atoms', natom), &
284        nctkdim_t('natom3', 3 * natom), &
285        nctkdim_t('number_of_phonon_modes', 3 * natom), &
286        nctkdim_t('anaddb_input_len', lenstr) &
287    ], defmode=.True.)
288    NCF_CHECK(ncerr)
289    ncerr = nctk_def_arrays(ana_ncid, [ &
290      nctkarr_t("anaddb_input_string", "char", "anaddb_input_len") &
291    ])
292    NCF_CHECK(ncerr)
293    !NCF_CHECK(nctk_defnwrite_ivars(ana_ncid, ["anaddb_version"], [1]))
294    NCF_CHECK(nctk_set_datamode(ana_ncid))
295    NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "anaddb_input_string"), string(:lenstr)))
296    NCF_CHECK(crystal_ncwrite(crystal, ana_ncid))
297 #endif
298  end if
299 
300  ! Calculation of Grunesein parameters.
301  if (inp%gruns_nddbs /= 0) then
302    call gruns_anaddb(inp, filnam(2), comm)
303    goto 50
304  end if
305 
306  ABI_ALLOCATE(instrain,(3*natom,6))
307  ABI_ALLOCATE(d2cart,(2,msize))
308  ABI_ALLOCATE(displ,(2*3*natom*3*natom))
309  ABI_ALLOCATE(eigval,(3,natom))
310  ABI_ALLOCATE(eigvec,(2,3,natom,3,natom))
311  ABI_ALLOCATE(phfrq,(3*natom))
312  ABI_ALLOCATE(zeff,(3,3,natom))
313  ABI_ALLOCATE(lst,(inp%nph2l))
314 
315 !**********************************************************************
316 !**********************************************************************
317 
318  ! Get Dielectric Tensor and Effective Charges
319  ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
320  iblok = ddb_get_dielt_zeff(ddb,crystal,inp%rfmeth,inp%chneut,inp%selectz,dielt,zeff)
321  ! Try to get dielt, in case just the DDE are present
322  if (iblok == 0) then
323    iblok_tmp = ddb_get_dielt(ddb,inp%rfmeth,dielt)
324  end if
325 
326  !if (iblok == 0) then
327  !  call wrtout(std_out, sjoin("- Cannot find dielectric tensor and Born effective charges in DDB file:", filnam(3)))
328  !  call wrtout(std_out, "Values initialized with zeros")
329  !else
330  !  call wrtout(std_out, sjoin("- Found dielectric tensor and Born effective charges in DDB file:", filnam(3)))
331  !end if
332 
333  if (my_rank == master) then
334 #ifdef HAVE_NETCDF
335    ! TODO: Cartesian or reduced?
336    ncerr = nctk_def_arrays(ana_ncid, [&
337    nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'),&
338    nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms")],&
339    defmode=.True.)
340    NCF_CHECK(ncerr)
341    ncerr = nctk_def_iscalars(ana_ncid, [character(len=nctk_slen) :: &
342        "asr", "chneut", "dipdip", "symdynmat"])
343    NCF_CHECK(ncerr)
344 
345    NCF_CHECK(nctk_set_datamode(ana_ncid))
346    NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, 'emacro_cart'), dielt))
347    NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, 'becs_cart'), zeff))
348    ncerr = nctk_write_iscalars(ana_ncid, [character(len=nctk_slen) :: &
349      "asr", "chneut", "dipdip", "symdynmat"], &
350      [inp%asr, inp%chneut, inp%dipdip, inp%symdynmat])
351    NCF_CHECK(ncerr)
352 #endif
353  end if
354 
355  ! Structural response at fixed polarization
356  if (inp%polflag == 1) then
357    ABI_ALLOCATE(d2flg,(msize))
358 
359    if(iblok/=0)then
360      ! Save the second-order derivatives
361      d2cart(1:2,1:msize) = ddb%val(1:2,1:msize,iblok)
362      d2flg(1:msize) = ddb%flg(1:msize,iblok)
363 
364    else
365      ! the gamma blok has not been found
366      if (inp%relaxat==0 .and. inp%relaxstr==0) then
367        ! The gamma blok is not needed
368        d2cart(1:2,1:msize)=zero
369        d2flg(1:msize)=1
370      else
371        ! There is a problem !
372        write(message, '(7a)' )&
373 &       'The dynamical matrix at Gamma is needed, in order to perform ',ch10,&
374 &       "relaxation at constant polarisation (Na Sai's method)",ch10,&
375 &       'However, this was not found in the DDB.',ch10,&
376 &       'Action: complete your DDB with the dynamical matrix at Gamma.'
377        MSG_ERROR(message)
378      end if
379    end if ! iblok not found
380 
381    ! Extract the block with the total energy
382    if (ddb_get_etotal(ddb,etotal) == 0) then
383      MSG_ERROR("DDB file does not contain GS etotal")
384    end if
385 
386    ! Extract the block with the gradients
387    ABI_ALLOCATE(fred,(3,natom))
388    qphon(:,:) = zero; qphnrm(:) = zero
389    rfphon(:) = 0; rfstrs(:) = 0; rfelfd(:) = 2
390    if (inp%relaxat == 1) rfphon(:) = 1
391    if (inp%relaxstr == 1) rfstrs(:) = 3
392 
393    call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp4)
394 
395    if (inp%relaxat == 1) then
396      index = 0
397      do iatom = 1, natom
398        do idir = 1, 3
399          index = index + 1
400          fred(idir,iatom) = ddb%val(1,index,iblok)
401        end do
402      end do
403    end if
404 
405    pel(1:3) = ddb%val(1,3*natom+4:3*natom+6,iblok)
406 
407    if (inp%relaxstr == 1) then
408      index = 3*natom + 6
409      do ii = 1, 6
410        index = index + 1
411        strten(ii) = ddb%val(1,index,iblok)
412      end do
413    end if
414 
415    ! when called from here red_ptot is not set! So set it to zero
416    red_ptot(:)=zero
417 
418    call relaxpol(Crystal,d2flg,d2cart,etotal,fred,inp%iatfix,&
419 &   ab_out,inp%istrfix,mpert,msize,inp%natfix,natom,&
420 &   inp%nstrfix,pel,red_ptot,inp%relaxat,inp%relaxstr,&
421 &   strten,inp%targetpol,usepaw)
422 
423    ABI_DEALLOCATE(fred)
424    ABI_DEALLOCATE(d2flg)
425  end if
426 
427 !***************************************************************************
428 
429  ! Compute non-linear optical susceptibilities and, if inp%nlflag < 3,
430  ! First-order change in the linear dielectric susceptibility induced by an atomic displacement
431  if (inp%nlflag > 0) then
432    ABI_ALLOCATE(dchide,(3,3,3))
433    ABI_ALLOCATE(dchidt,(natom,3,3,3))
434 
435    if (ddb_get_dchidet(ddb,inp%ramansr,inp%nlflag,dchide,dchidt) == 0) then
436      MSG_ERROR("Cannot find block corresponding to non-linear optical susceptibilities in DDB file")
437    end if
438 
439    ! Save to the netcdf
440    if (my_rank == master) then
441 #ifdef HAVE_NETCDF
442      ncerr = nctk_def_arrays(ana_ncid, [nctkarr_t("dchide", "dp", "three, three, three")], defmode=.True.)
443      NCF_CHECK(ncerr)
444      NCF_CHECK(nctk_set_datamode(ana_ncid))
445      NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "dchide"), dchide))
446 
447      ! dchidt only present if nlflag==1 or 2
448      if (inp%nlflag < 3) then
449        ncerr = nctk_def_arrays(ana_ncid, [nctkarr_t("dchidt", "dp", &
450        "number_of_atoms, three, three, three")], defmode=.True.)
451        NCF_CHECK(ncerr)
452        NCF_CHECK(nctk_set_datamode(ana_ncid))
453        NCF_CHECK(nf90_put_var(ana_ncid, nctk_idname(ana_ncid, "dchidt"), dchidt))
454      end if
455 #endif
456    end if
457 
458  end if ! nlflag
459 
460 !**********************************************************************
461 ! Interatomic Forces Calculation
462 !**********************************************************************
463  if (inp%ifcflag ==1) then
464    ! ifc to be calculated for interpolation
465    write(message, '(a,a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,&
466 &   ' Calculation of the interatomic forces ',ch10
467    call wrtout(std_out,message,'COLL')
468    call wrtout(ab_out,message,'COLL')
469 
470 ! TODO : check if this wrtout should be removed in latest merge 17 feb 2017
471    call timein(tcpu,twall)
472    write(message, '(a,f11.3,a,f11.3,a)' )'-begin at tcpu',tcpu-tcpui,'  and twall',twall-twalli,' sec'
473    call wrtout(std_out,message,'COLL')
474    call wrtout(ab_out,message,'COLL')
475 
476    if (any(inp%qrefine(:) > 1)) then
477      ! Gaal-Nagy's algorithm in PRB 73 014117 [[cite:GaalNagy2006]]
478 
479      ! Build the IFCs using the coarse q-mesh.
480      do ii = 1, 3
481        ngqpt_coarse(ii) = inp%ngqpt(ii)/inp%qrefine(ii)
482      end do
483      call ifc_init(Ifc_coarse,Crystal,ddb,&
484 &     inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,ngqpt_coarse,inp%nqshft,inp%q1shft,dielt,zeff,&
485 &     inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm)
486 
487      ! And now use the coarse q-mesh to fill the entries in dynmat(q)
488      ! on the dense q-mesh that cannot be obtained from the DDB file.
489      call ifc_init(Ifc,Crystal,ddb,&
490 &     inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,zeff,&
491 &     inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm,Ifc_coarse=Ifc_coarse)
492      call ifc_free(Ifc_coarse)
493 
494    else
495      call ifc_init(Ifc,Crystal,ddb,&
496 &     inp%brav,inp%asr,inp%symdynmat,inp%dipdip,inp%rfmeth,inp%ngqpt(1:3),inp%nqshft,inp%q1shft,dielt,zeff,&
497 &     inp%nsphere,inp%rifcsph,inp%prtsrlr,inp%enunit,comm)
498    end if
499 
500    call ifc_print(ifc, unit=std_out)
501 
502    ! Compute speed of sound.
503    if (inp%vs_qrad_tolkms(1) > zero) then
504      call ifc_speedofsound(ifc, crystal, inp%vs_qrad_tolkms, ana_ncid, comm)
505      !call ifc_test_phinterp(ifc, crystal, [8,8,8], 1, [zero,zero,zero], [3,3,3], comm, test_dwdq=.True.)
506      !stop
507    end if
508 
509    ! Print analysis of the real-space interatomic force constants
510    ! TODO: ifc_out should not have side effects
511    if (my_rank == master .and. inp%ifcout/=0) then
512      call ifc_write(Ifc,inp%ifcana,inp%atifc,inp%ifcout,inp%prt_ifc,ana_ncid)
513    end if
514  end if
515 
516 !**********************************************************************
517 
518 !Electron-phonon section
519  if (inp%elphflag == 1) then
520    call elphon(inp,Crystal,Ifc,filnam,comm)
521  end if
522 
523 !**********************************************************************
524 
525  if (sum(abs(inp%thermal_supercell))>0 .and. inp%ifcflag==1) then
526    ABI_ALLOCATE(thm_scells, (inp%ntemper))
527    call zacharias_supercell_make(Crystal, Ifc, inp%ntemper, inp%thermal_supercell, inp%tempermin, inp%temperinc, thm_scells)
528    call zacharias_supercell_print(filnam(2), inp%ntemper, inp%tempermin, inp%temperinc, thm_scells)
529  end if
530 
531  ! Phonon density of states calculation, Start if interatomic forces have been calculated
532  if (inp%ifcflag==1 .and. any(inp%prtdos==[1, 2])) then
533    write(message,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,' Calculation of phonon density of states ',ch10
534    call wrtout(ab_out,message,'COLL')
535    call wrtout(std_out,message,'COLL')
536 
537    ! Only 1 shift in q-mesh
538    wminmax = zero
539    do
540      call mkphdos(Phdos, Crystal, Ifc, inp%prtdos, inp%dosdeltae, inp%dossmear, inp%ng2qpt, 1, inp%q2shft, &
541      wminmax, count_wminmax, comm)
542      if (all(count_wminmax == 0)) exit
543      wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05
544      wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05
545      call phdos_free(phdos)
546    end do
547 
548    if (iam_master) then
549      call phdos_print_msqd(Phdos, filnam(2), inp%ntemper, inp%tempermin, inp%temperinc)
550      call phdos_print(Phdos, strcat(filnam(2), "_PHDOS"))
551      call phdos_print_debye(Phdos, Crystal%ucvol)
552      call phdos_print_thermo(PHdos, strcat(filnam(2), "_THERMO"), inp%ntemper, inp%tempermin, inp%temperinc)
553 
554 #ifdef HAVE_NETCDF
555      ncerr = nctk_open_create(phdos_ncid, strcat(filnam(2), "_PHDOS.nc"), xmpi_comm_self)
556      NCF_CHECK_MSG(ncerr, "Creating PHDOS.nc file")
557      NCF_CHECK(crystal_ncwrite(Crystal, phdos_ncid))
558      call phdos_ncwrite(Phdos, phdos_ncid)
559      NCF_CHECK(nf90_close(phdos_ncid))
560 #endif
561    end if
562 
563    call phdos_free(Phdos)
564  end if
565 
566  if (iam_master .and. inp%ifcflag==1 .and. inp%outboltztrap==1) then
567    call ifc_outphbtrap(Ifc,Crystal,inp%ng2qpt,1,inp%q2shft,filnam(2))
568  end if
569 
570  ! Phonon density of states and thermodynamical properties calculation
571  ! Start if interatomic forces and thermal flags are on
572  if (inp%ifcflag==1 .and. any(inp%thmflag==[1,2])) then
573 
574    write(message, '(a,(80a),a,a,a,a,a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,&
575 &   ' Calculation of phonon density of states, ',ch10,&
576 &   '    thermodynamical properties, ',ch10,&
577 &   '    and Debye-Waller factors.',ch10
578    call wrtout(ab_out,message,'COLL')
579    call wrtout(std_out,message,'COLL')
580 
581    if (inp%thmflag==1) then
582      call harmonic_thermo(Ifc,Crystal,ddb%amu,inp,ab_out,filnam(2),comm)
583 
584    else if (inp%thmflag==2) then
585      write(message, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Entering thm9 routine with thmflag=2 ',ch10
586      call wrtout(std_out,message,'COLL')
587      call wrtout(ab_out,message,'COLL')
588 
589      call harmonic_thermo(Ifc,Crystal,ddb%amu,inp,ab_out,filnam(2),comm,thmflag=inp%thmflag)
590    end if
591  end if
592 
593 !**********************************************************************
594 
595  ! Now treat the first list of vectors (without non-analyticities)
596  call mkphbs(Ifc,Crystal,inp,ddb,asrq0,filnam(2),comm)
597 
598 !***********************************************************************
599 
600  ! Interpolate the DDB onto the first list of vectors and write the file.
601 
602  if (inp%prtddb==1 .and. inp%ifcflag==1) then
603 
604    call ddb_hdr_open_read(ddb_hdr,filnam(3),ddbun,DDB_VERSION)
605    close(ddbun)
606 
607    call ddb_interpolate(Ifc,Crystal,inp,ddb,ddb_hdr,asrq0,filnam(2),comm)
608 
609    call ddb_hdr_free(ddb_hdr)
610  end if
611 
612 !***********************************************************************
613 
614  if (inp%thmflag>=3 .and. inp%thmflag<=8) then
615    !write(std_out,*)'Entering thmeig: '
616    elph_base_name=trim(filnam(2))//"_ep"
617 
618    call thmeig(inp,ddb,Crystal,elph_base_name,filnam(5),&
619 &   ddbun,ab_out,natom,mpert,msize,asrq0%d2asr,comm)
620  end if
621 
622 !**********************************************************************
623 
624 !Now treat the second list of vectors (only at the Gamma point,
625 !but can include non-analyticities), as well as the frequency-dependent dielectric tensor
626 
627  if (inp%nlflag > 0) then
628    ABI_ALLOCATE(rsus,(3*natom,3,3))
629  end if
630  ABI_ALLOCATE(fact_oscstr,(2,3,3*natom))
631 
632  if (inp%nph2l/=0 .or. inp%dieflag==1) then
633 
634    write(message, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,&
635 &   ch10,' Treat the second list of vectors ',ch10
636    call wrtout(std_out,message,'COLL')
637    call wrtout(ab_out,message,'COLL')
638 
639    ! Before examining every direction or the dielectric tensor, generates the dynamical matrix at gamma
640    qphon(:,1)=zero; qphnrm(1)=zero
641 
642    ! Generation of the dynamical matrix in cartesian coordinates
643    if (inp%ifcflag==1) then
644 
645      ! Get d2cart using the interatomic forces and the
646      ! long-range coulomb interaction through Ewald summation
647      call gtdyn9(ddb%acell,Ifc%atmfrc,dielt,inp%dipdip,&
648 &     Ifc%dyewq0,d2cart,Crystal%gmet,ddb%gprim,mpert,natom,&
649 &     Ifc%nrpt,qphnrm(1),qphon,Crystal%rmet,ddb%rprim,Ifc%rpt,&
650 &     Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,zeff)
651 
652    else if (inp%ifcflag==0) then
653 
654      ! Look after the information in the DDB
655      rfphon(1:2)=1; rfelfd(1:2)=2; rfstrs(1:2)=0
656      call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,inp%rfmeth)
657 
658      ! Copy the dynamical matrix in d2cart
659      d2cart(:,1:msize)=ddb%val(:,:,iblok)
660 
661      ! Eventually impose the acoustic sum rule
662      call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
663    end if ! end of the generation of the dynamical matrix at gamma.
664 
665    if (inp%nph2l/=0) then
666 
667      if (my_rank == master) then
668 #ifdef HAVE_NETCDF
669        iphl2 = 0
670        call nctk_defwrite_nonana_terms(ana_ncid, iphl2, inp%nph2l, inp%qph2l, natom, phfrq, displ, "define")
671 #endif
672      end if
673 
674      ! Examine every wavevector of this list
675      do iphl2=1,inp%nph2l
676 
677        ! Initialisation of the phonon wavevector
678        qphon(:,1)=inp%qph2l(:,iphl2)
679        qphnrm(1)=inp%qnrml2(iphl2)
680 
681        ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
682        call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,&
683 &       mpert,msym,natom,nsym,ntypat,phfrq,qphnrm(1),qphon,Crystal%rprimd,inp%symdynmat,&
684 &       Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
685 
686        ! Write the phonon frequencies
687        call dfpt_prtph(displ,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon)
688 
689        if (my_rank == master) then
690 #ifdef HAVE_NETCDF
691          ! Loop is not MPI-parallelized --> no need for MPI-IO API.
692          call nctk_defwrite_nonana_terms(ana_ncid, iphl2, inp%nph2l, inp%qph2l, natom, phfrq, displ, "write")
693 #endif
694        end if
695 
696        ! Determine the symmetries of the phonon modes at Gamma
697        if (sum(abs(qphon(:,1)))<DDB_QTOL) then
698          call dfpt_symph(ab_out,ddb%acell,eigvec,Crystal%indsym,natom,nsym,phfrq,ddb%rprim,Crystal%symrel)
699        end if
700 
701        ! Write Raman susceptibilities
702        if (inp%nlflag == 1) then
703          call ramansus(d2cart,dchide,dchidt,displ,mpert,natom,phfrq,qphon,qphnrm(1),rsus,Crystal%ucvol)
704        end if
705 
706        ! Prepare the evaluation of the Lyddane-Sachs-Teller relation
707        if(inp%dieflag==1 .and. natom>1)then
708          lst(iphl2)=zero
709          ! The fourth mode should have positive frequency, otherwise,
710          ! there is an instability, and the LST relationship should not be evaluated
711          if(phfrq(4)>tol6)then
712            do ii=4,3*natom
713              lst(iphl2)=lst(iphl2)+2*log(phfrq(ii))
714            end do
715          end if
716        end if
717 
718      end do ! iphl2
719    end if ! nph2l/=0
720 
721    ! The frequency-dependent dielectric tensor (and oscillator strength).
722    if (inp%dieflag==1)then
723      write(message, '(6a)' )&
724 &     ' the frequency-dependent dielectric tensor (and also once more',ch10,&
725 &     ' the phonons at gamma - without non-analytic part )',ch10,ch10,&
726 &     ' The frequency-dependent dielectric tensor'
727      call wrtout(std_out,message,'COLL')
728 
729      ! Initialisation of the phonon wavevector
730      qphon(:,1)=zero; qphnrm(1)=zero
731 
732      ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix
733      call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,&
734 &     mpert,msym,natom,nsym,ntypat,phfrq,qphnrm(1),qphon,&
735 &     Crystal%rprimd,inp%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
736 
737      ! Write the phonon frequencies (not to ab_out, however)
738      call dfpt_prtph(displ,0,inp%enunit,-1,natom,phfrq,qphnrm(1),qphon)
739 
740      ! Evaluation of the oscillator strengths and frequency-dependent dielectric tensor.
741      call ddb_diel(Crystal,ddb%amu,inp,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
742 &     ab_out,lst,mpert,natom,inp%nph2l,phfrq,comm,ana_ncid)
743      ! write(std_out,*)'after ddb_diel, dielt_rlx(:,:)=',dielt_rlx(:,:)
744    end if
745 
746    ! If the electronic dielectric tensor only is needed...
747    if (inp%dieflag==2.or.inp%dieflag==3.or. inp%dieflag==4) then
748 !    Everything is already in place...
749      call ddb_diel(Crystal,ddb%amu,inp,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
750 &     ab_out,lst,mpert,natom,inp%nph2l,phfrq,comm,ana_ncid)
751    end if
752 
753  end if ! either nph2l/=0  or  dieflag==1
754 
755 !**********************************************************************
756 
757 !In case nph2l was equal to 0, the electronic dielectric tensor has to be computed independently.
758 
759  if (inp%dieflag==2 .and. inp%nph2l==0) then
760    call wrtout(std_out,'nph2l=0, so compute the electronic dielectric tensor independently','COLL')
761    ! Look after the second derivative matrix at gamma in the DDB
762    ! Note that the information on the dielectric tensor is completely
763    ! independent of the interatomic force constant calculation
764    qphon(:,1)=zero; qphnrm(1)=zero
765    rfphon(1:2)=0; rfelfd(1:2)=2; rfstrs(1:2)=0
766    call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,inp%rfmeth)
767    if (iblok == 0) then
768      MSG_ERROR("DDB file must contain both derivatives wrt electric field, Check your calculations")
769    end if
770 
771    d2cart(:,1:msize)=ddb%val(:,:,iblok)
772 
773    ! Print the electronic dielectric tensor
774    call ddb_diel(Crystal,ddb%amu,inp,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
775    ab_out,lst,mpert,natom,inp%nph2l,phfrq,comm,ana_ncid)
776  end if
777 
778 !**********************************************************************
779 
780  ! Compute the electrooptic tensor
781  if (inp%nlflag == 1) then
782    ! In case dieflag = 2, recompute phonon frequencies and eigenvectors without non-analyticity
783    if (inp%dieflag == 2) then
784      qphon(:,1)=zero; qphnrm(1)=zero
785      call dfpt_phfrq(ddb%amu,displ,d2cart,eigval,eigvec,Crystal%indsym,&
786 &     mpert,msym,natom,nsym,ntypat,phfrq,qphnrm(1),qphon,&
787 &     Crystal%rprimd,inp%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol)
788    end if
789 
790    rsus = zero
791    call ramansus(d2cart,dchide,dchidt,displ,mpert,natom,phfrq(1),qphon,qphnrm(1),rsus,Crystal%ucvol)
792 
793    call electrooptic(dchide,inp%dieflag,epsinf,fact_oscstr,natom,phfrq,inp%prtmbm,rsus,Crystal%ucvol)
794  end if ! condition on nlflag
795 
796  ABI_DEALLOCATE(fact_oscstr)
797  if (inp%nlflag > 0) then
798    ABI_DEALLOCATE(dchide)
799    ABI_DEALLOCATE(rsus)
800    ABI_DEALLOCATE(dchidt)
801  end if
802 
803 !**********************************************************************
804 
805 ! Here treating the internal strain tensors at Gamma point
806  if (inp%instrflag/=0) then
807 
808    write(message, '(a,a,(80a),a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
809 &   ' Calculation of the internal-strain  tensor',ch10
810    call wrtout(std_out,message,'COLL')
811    call wrtout(ab_out,message,'COLL')
812 
813    if (inp%instrflag==1) then
814      call wrtout(std_out,'instrflag=1, so extract the internal strain constant from the 2DTE','COLL')
815 
816      ! looking after the no. of blok that contains the internal strain tensor
817      qphon(:,1)=zero; qphnrm(1)=zero
818      rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
819 
820      call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,inp%rfmeth)
821      if (iblok == 0) then
822        MSG_ERROR("DDB file must contain both uniaxial and shear strain for piezoelectric, Check your calculations")
823      end if
824 
825      ! then print the internal stain tensor
826      prt_internalstr=2
827      call ddb_internalstr(inp%asr,ddb%val,asrq0%d2asr,iblok,instrain,ab_out,mpert,natom,ddb%nblok,prt_internalstr)
828    end if
829  end if !end the part for internal strain
830 
831 !**********************************************************************
832 
833 !here treating the elastic tensors at Gamma Point
834  if (inp%elaflag/=0) then
835    write(message, '(a,a,(80a),a,a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
836 &   ' Calculation of the elastic and compliances tensor (Voigt notation)',ch10
837    call wrtout(std_out,message,'COLL')
838    call wrtout(ab_out,message,'COLL')
839 
840    if (any(inp%elaflag == [1,2,3,4,5])) then
841      call wrtout(std_out,'so extract the elastic constant from the 2DTE','COLL')
842 
843      ! look after the blok no. that contains the stress tensor
844      qphon(:,1)=zero; qphnrm(1)=zero
845      rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=0
846 
847      call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp4)
848      iblok_stress=iblok
849 
850      ! look after the blok no.iblok that contains the elastic tensor
851      qphon(:,1)=zero; qphnrm(1)=zero
852      rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
853 
854      ! for both diagonal and shear parts
855      call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,inp%rfmeth)
856      if (iblok == 0) then
857        MSG_ERROR("DDB file must contain both uniaxial and shear strain when elaflag != 0, Check your calculations")
858      end if
859 
860      ! print the elastic tensor
861      call ddb_elast(inp,crystal,ddb%val,compl,compl_clamped,compl_stress,asrq0%d2asr,&
862 &     elast,elast_clamped,elast_stress,iblok,iblok_stress,&
863 &     instrain,ab_out,mpert,natom,ddb%nblok,ana_ncid)
864    end if
865  end if !ending the part for elastic tensors
866 
867 !**********************************************************************
868 
869 !here treating the piezoelectric tensor at Gamma Point
870  if (inp%piezoflag/=0 .or. inp%dieflag==4 .or. inp%elaflag==4) then
871    write(message, '(a,a,(80a),a,a,a,a,a)') ch10,('=',ii=1,80),ch10,ch10,&
872 &   ' Calculation of the tensor related to piezoelectric effetc',ch10,&
873 &   '  (Elastic indices in Voigt notation)',ch10
874    call wrtout(std_out,message,'COLL')
875    call wrtout(ab_out,message,'COLL')
876 
877    if (any(inp%piezoflag == [1,2,3,4,5,6,7]) .or. inp%dieflag==4 .or.inp%elaflag==4) then
878      call wrtout(std_out,'extract the piezoelectric constant from the 2DTE','COLL')
879 
880      ! looking for the gamma point block
881      qphon(:,1)=zero; qphnrm(1)=zero
882      rfphon(1:2)=0; rfelfd(1:2)=0; rfstrs(1:2)=3
883 
884      ! for both diagonal and shear parts
885      call gtblk9(ddb,iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,inp%rfmeth)
886      if (iblok == 0) then
887        MSG_ERROR("DDB file must contain both uniaxial and shear strain for piezoelectric, Check your calculations")
888      end if
889 
890      ! then print out the piezoelectric constants
891      call ddb_piezo(inp,ddb%val,dielt_rlx,elast,iblok,instrain,ab_out,mpert,natom,ddb%nblok,piezo,Crystal%ucvol,ana_ncid)
892    end if
893  end if
894 
895 !**********************************************************************
896 
897  ! Free memory
898  ABI_DEALLOCATE(displ)
899  ABI_DEALLOCATE(d2cart)
900  ABI_DEALLOCATE(eigval)
901  ABI_DEALLOCATE(eigvec)
902  ABI_DEALLOCATE(lst)
903  ABI_DEALLOCATE(phfrq)
904  ABI_DEALLOCATE(zeff)
905  ABI_DEALLOCATE(instrain)
906 
907  50 continue
908 
909  call asrq0_free(asrq0)
910  call ifc_free(Ifc)
911  call crystal_free(Crystal)
912  call ddb_free(ddb)
913  call anaddb_dtset_free(inp)
914  call thermal_supercell_free(inp%ntemper, thm_scells)
915 
916  if (sum(abs(inp%thermal_supercell))>0 .and. inp%ifcflag==1) then
917    ABI_DEALLOCATE(thm_scells)
918  end if
919 
920  ! Close files
921  if (iam_master) then
922 #ifdef HAVE_NETCDF
923    NCF_CHECK(nf90_close(ana_ncid))
924 #endif
925  end if
926 
927  call timein(tcpu,twall)
928  tsec(1)=tcpu-tcpui; tsec(2)=twall-twalli
929  write(message, '(a,i4,a,f13.1,a,f13.1)' )' Proc.',my_rank,' individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
930  call wrtout(std_out,message,"COLL")
931 
932  if (iam_master) then
933    write(ab_out, '(a,a,a,i4,a,f13.1,a,f13.1)' )'-',ch10,&
934 &   '- Proc.',my_rank,' individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
935  end if
936 
937  call xmpi_sum(tsec,comm,ierr)
938 
939  write(message, '(a,(80a),a,a,a,f11.3,a,f11.3,a,a,a,a)' ) ch10,&
940 & ('=',ii=1,80),ch10,ch10,&
941 & '+Total cpu time',tsec(1),&
942 & '  and wall time',tsec(2),' sec',ch10,ch10,&
943 & ' anaddb : the run completed succesfully.'
944  call wrtout(std_out,message,'COLL')
945  call wrtout(ab_out,message,'COLL')
946 
947  if (iam_master) then
948    ! Write YAML document with the final summary.
949    ! we use this doc to test whether the calculation is completed.
950    write(std_out,"(a)")"--- !FinalSummary"
951    write(std_out,"(a)")"program: anaddb"
952    write(std_out,"(2a)")"version: ",trim(abinit_version)
953    write(std_out,"(2a)")"start_datetime: ",start_datetime
954    write(std_out,"(2a)")"end_datetime: ",asctime()
955    write(std_out,"(a,f13.1)")"overall_cpu_time: ",tsec(1)
956    write(std_out,"(a,f13.1)")"overall_wall_time: ",tsec(2)
957    write(std_out,"(a,i0)")"mpi_procs: ",xmpi_comm_size(xmpi_world)
958    write(std_out,"(a,i0)")"omp_threads: ",xomp_get_num_threads(open_parallel=.True.)
959    !write(std_out,"(a,i0)")"num_warnings: ",nwarning
960    !write(std_out,"(a,i0)")"num_comments: ",ncomment
961    write(std_out,"(a)")"..."
962    call flush_unit(std_out)
963  end if
964 
965  ! Write information on file about the memory before ending mpi module, if memory profiling is enabled
966  call abinit_doctor(filnam(2))
967 
968  call flush_unit(ab_out)
969  call flush_unit(std_out)
970 
971  if (iam_master) close(ab_out)
972 
973  100 call xmpi_end()
974 
975  end program anaddb