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