TABLE OF CONTENTS
ABINIT/atdep [ Programs ]
NAME
atdep
FUNCTION
Calculations of phonons using molecular dynamic simulations
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (FB,JB) 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
NOTES
The input files are input.in, xred.dat, fcart.dat and etot.dat See the examples in the test directory
TODO
A lot of things
INPUTS
(main routine)
OUTPUT
(main routine)
SOURCE
31 #if defined HAVE_CONFIG_H 32 #include "config.h" 33 #endif 34 35 #include "abi_common.h" 36 37 program atdep 38 39 use defs_basis 40 use m_abicore 41 use m_phonons 42 use m_errors 43 use m_abi_linalg 44 use m_xmpi 45 use m_abihist 46 use m_io_tools 47 use m_argparse 48 49 use m_ifc, only : ifc_type 50 use m_crystal, only : crystal_t, crystal_free 51 use m_ddb, only : ddb_type 52 use m_tdep_abitypes, only : Qbz_type, tdep_init_crystal, tdep_init_ifc, tdep_init_ddb, tdep_write_ddb, & 53 & tdep_destroy_qbz, tdep_destroy_ddb, tdep_ifc2phi2 54 use m_tdep_phi4, only : tdep_calc_phi4fcoeff, tdep_calc_phi4ref, tdep_write_phi4, tdep_calc_ftot4 55 use m_tdep_phi3, only : tdep_calc_phi3fcoeff, tdep_calc_phi3ref, tdep_write_phi3, tdep_calc_ftot3, & 56 & tdep_calc_alpha_gamma, tdep_write_gruneisen 57 use m_tdep_phi2, only : tdep_calc_phi2fcoeff, tdep_calc_phi1fcoeff, tdep_calc_phi2, tdep_write_phi2, tdep_calc_ftot2, & 58 & Eigen_type, tdep_init_eigen2nd, tdep_destroy_eigen2nd, tdep_calc_phi1, & 59 & tdep_write_phi1, tdep_init_phi2, tdep_destroy_phi2, Phi2_type 60 use m_tdep_latt, only : tdep_make_latt, Lattice_type 61 use m_tdep_sym, only : tdep_make_sym, Symetries_type, tdep_destroy_sym 62 use m_tdep_readwrite, only : tdep_print_Aknowledgments, tdep_read_input, tdep_distrib_data, tdep_init_MPIdata, & 63 & tdep_destroy_mpidata, Input_type, MPI_enreg_type, tdep_destroy_invar 64 use m_tdep_utils, only : Coeff_Moore_type, tdep_calc_MoorePenrose, tdep_MatchIdeal2Average, tdep_calc_model 65 use m_tdep_qpt, only : tdep_make_qptpath, Qpoints_type, tdep_destroy_qpt 66 use m_tdep_phdos, only : tdep_calc_phdos,tdep_calc_elastic,tdep_calc_thermo 67 use m_tdep_shell, only : Shell_type, tdep_init_shell2at, tdep_init_shell3at, tdep_init_shell4at, & 68 & tdep_init_shell1at, tdep_destroy_shell 69 use m_tdep_constraints, only : tdep_calc_constraints, tdep_check_constraints 70 71 implicit none 72 73 integer :: natom,natom_unitcell,ncoeff1st,ncoeff2nd,ncoeff3rd,ncoeff4th,ntotcoeff,ntotconst 74 integer :: stdout,stdlog,nshell_max,ii,jj,ishell,istep,iatom 75 integer :: print_mem_report,ierr 76 double precision :: U0 77 double precision, allocatable :: ucart(:,:,:),proj1st(:,:,:),proj2nd(:,:,:),proj3rd(:,:,:),proj4th(:,:,:) 78 double precision, allocatable :: proj_tmp(:,:,:),Forces_TDEP(:),Fresid(:) 79 double precision, allocatable :: Phi1(:) ,Phi1_coeff(:,:) 80 double precision, allocatable :: Phi2_coeff(:,:) 81 double precision, allocatable :: Phi3_ref(:,:,:,:) ,Phi3_coeff(:,:) 82 double precision, allocatable :: Phi4_ref(:,:,:,:,:),Phi4_coeff(:,:) 83 double precision, allocatable :: Forces_MD(:),MP_coeff(:,:) 84 double precision, allocatable :: distance(:,:,:),Rlatt_cart(:,:,:),Rlatt4Abi(:,:,:) 85 double precision, allocatable :: Phi1Ui(:),Phi2UiUj(:),Phi3UiUjUk(:),Phi4UiUjUkUl(:) 86 type(args_t) :: args 87 type(phdos_t) :: PHdos 88 type(Phi2_type) :: Phi2 89 type(Input_type) :: Invar 90 type(Lattice_type) :: Lattice 91 type(Symetries_type) :: Sym 92 type(Qpoints_type) :: Qpt 93 type(Qbz_type) :: Qbz 94 type(ifc_type) :: Ifc 95 type(ddb_type) :: DDB 96 type(crystal_t) :: Crystal 97 type(Shell_type) :: Shell1at,Shell2at,Shell3at,Shell4at 98 type(Coeff_Moore_type) :: CoeffMoore 99 type(Eigen_type) :: Eigen2nd_MP 100 type(Eigen_type) :: Eigen2nd_path 101 type(MPI_enreg_type) :: MPIdata 102 type(abihist) :: Hist 103 104 !TEST 105 !FB integer, allocatable :: data_tmp(:,:),data_loc(:,:),data_gather(:,:),shft_step(:) 106 !FB integer, allocatable :: nstep_acc(:),tab_step(:) 107 !FB integer :: istep,ierr,iproc,dim1,dim2,remain,idim1 108 !TEST 109 110 !****************************************************************** 111 112 !========================================================================================== 113 !===================== Initialization & Reading ========================================== 114 !========================================================================================== 115 ! Change communicator for I/O (mandatory!) 116 call abi_io_redirect(new_io_comm=xmpi_world) 117 ! Initialize MPI 118 call xmpi_init() 119 120 ! Parse command line arguments. 121 args = args_parser(); if (args%exit /= 0) goto 100 122 123 ! Initialize memory profiling if activated at configure time. 124 ! if a full report is desired, set the argument of abimem_init to "2" instead of "0" via the command line. 125 ! note that the file can easily be multiple GB in size so don't use this option normally 126 #ifdef HAVE_MEM_PROFILING 127 call abimem_init(args%abimem_level, limit_mb=args%abimem_limit_mb) 128 #endif 129 130 ! Read input values from the input.in input file 131 call tdep_read_input(Hist,Invar) 132 call tdep_init_MPIdata(Invar,MPIdata) 133 call tdep_distrib_data(Hist,Invar,MPIdata) 134 call abihist_free(Hist) 135 136 !FB TEST 137 !FB Invar%nstep_tot=10 138 !FB dim1=1 139 !FB dim2=1 140 !FB ABI_MALLOC(data_tmp,(dim1,Invar%nstep_tot*dim2)); data_tmp(:,:)=zero 141 !FB do istep=1,Invar%nstep_tot 142 !FB data_tmp(:,dim2*(istep-1)+1:dim1*istep)=istep 143 !FB end do 144 !FB if (MPIdata%iam_master) write(Invar%stdlog,*) MPIdata%me_step,data_tmp(:,:) 145 !FB 146 !FB Invar%my_nstep=int(Invar%nstep_tot/MPIdata%nproc) 147 !FB remain=Invar%nstep_tot-MPIdata%nproc*Invar%my_nstep 148 !FB do ii=1,remain 149 !FB if ((ii-1).eq.MPIdata%me_step) Invar%my_nstep=Invar%my_nstep+1 150 !FB end do 151 !FB!FB ABI_MALLOC(MPIdata%nstep_all,(MPIdata%nproc)); MPIdata%nstep_all(:)=zero 152 !FB call xmpi_allgather(Invar%my_nstep,MPIdata%nstep_all,MPIdata%comm_step,ierr) 153 !FB if (MPIdata%me_step.eq.MPIdata%master) write(Invar%stdlog,*) 'Nstep(nproc)=',MPIdata%nstep_all(:) 154 !FB call flush_unit(6) 155 !FB call xmpi_barrier(MPIdata%comm_step) 156 !FB 157 !FB ABI_MALLOC(nstep_acc,(MPIdata%nproc+1)); nstep_acc(:)=zero 158 !FB nstep_acc(1)=0 159 !FB do ii=2,MPIdata%nproc+1 160 !FB nstep_acc(ii)=nstep_acc(ii-1)+MPIdata%nstep_all(ii-1) 161 !FB end do 162 !FB if (MPIdata%me_step.eq.MPIdata%master) write(Invar%stdlog,*) 'NSTEP_ACC=',nstep_acc(:) 163 !FB call flush_unit(6) 164 !FB call xmpi_barrier(MPIdata%comm_step) 165 !FB 166 !FB ABI_MALLOC(tab_step,(Invar%nstep_tot)); tab_step(:)=zero 167 !FB! First distrib 168 !FB do iproc=1,MPIdata%nproc 169 !FB do istep=1,Invar%nstep_tot 170 !FB if ((istep.gt.nstep_acc(iproc)).and.(istep.le.nstep_acc(iproc+1))) then 171 !FB tab_step(istep)=iproc-1 172 !FB end if 173 !FB end do 174 !FB end do 175 !FB 176 !FB! Second distrib 177 !FB tab_step(:)=zero 178 !FB do istep=1,Invar%nstep_tot 179 !FB do iproc=1,MPIdata%nproc 180 !FB if (mod((istep-1),MPIdata%nproc).eq.(iproc-1)) then 181 !FB tab_step(istep)=iproc-1 182 !FB end if 183 !FB end do 184 !FB end do 185 !FB 186 !FB if (MPIdata%me_step.eq.MPIdata%master) write(Invar%stdlog,*) 'TAB_STEP=',tab_step(:) 187 !FB call flush_unit(6) 188 !FB call xmpi_barrier(MPIdata%comm_step) 189 !FB 190 !FB if (nstep_acc(MPIdata%nproc+1).ne.Invar%nstep_tot) then 191 !FB write(Invar%stdlog,*) 'STOP : pb in nstep_acc' 192 !FB stop 193 !FB end if 194 !FB 195 !FB ii=1 196 !FB ABI_MALLOC(data_loc,(dim1,Invar%my_nstep*dim2)); data_loc(:,:)=zero 197 !FB do istep=1,Invar%nstep_tot 198 !FB if (tab_step(istep).eq.MPIdata%me_step) then 199 !FB data_loc(:,dim2*(ii-1)+1:dim2*ii)=data_tmp(:,dim2*(istep-1)+1:dim2*istep) 200 !FB ii=ii+1 201 !FB end if 202 !FB end do 203 !FB 204 !FB do iproc=1,MPIdata%nproc 205 !FB if (iproc-1.eq.MPIdata%me_step) write(Invar%stdlog,*) 'DATA_LOC =',MPIdata%me_step,data_loc(:,:) 206 !FB end do 207 !FB call flush_unit(6) 208 !FB call xmpi_barrier(MPIdata%comm_step) 209 !FB 210 !FB!FB ABI_MALLOC(shft_step,(MPIdata%nproc)); shft_step(:)=zero 211 !FB!FB shft_step(1)=0 212 !FB!FB do ii=2,MPIdata%nproc 213 !FB!FB shft_step(ii)=shft_step(ii-1)+MPIdata%nstep_all(ii-1) 214 !FB!FB end do 215 !FB!FB 216 !FB!FB ABI_MALLOC(data_gather,(dim1,Invar%nstep_tot*dim2)); data_gather(:,:)=zero 217 !FB!FB do idim1=1,dim1 218 !FB!FB call xmpi_gatherv(data_loc(idim1,:),Invar%my_nstep*dim2,data_gather(idim1,:),MPIdata%nstep_all*dim2,dim2*shft_step,& 219 !FB!FB& MPIdata%master,MPIdata%comm_step,ierr) 220 !FB!FB end do 221 !FB!FB 222 !FB!FB if (MPIdata%me_step.eq.MPIdata%master) write(Invar%stdlog,*) "============================================" 223 !FB!FB if (MPIdata%me_step.eq.MPIdata%master) write(Invar%stdlog,*) MPIdata%me_step,data_gather(:,:) 224 !FB!FB call flush_unit(6) 225 !FB!!FBFB 226 !FB 227 !FB ABI_MALLOC(data_gather,(dim1,Invar%nstep_tot*dim2)); data_gather(:,:)=zero 228 !FB call xmpi_allgatherv(data_loc,dim1*Invar%my_nstep,data_gather,dim1*MPIdata%nstep_all,dim1*nstep_acc(1:MPIdata%nproc),MPIdata%comm_step,ierr) 229 !FB!FB call xmpi_scatterv(data_gather,MPIdata%nstep_all,nstep_acc,data_loc,Invar%my_nstep,& 230 !FB!FB& MPIdata%master,MPIdata%comm_step,ierr) 231 !FB do iproc=1,MPIdata%nproc 232 !FB if (MPIdata%me_step.eq.(iproc-1)) write(Invar%stdlog,*) "============================================" 233 !FB if (MPIdata%me_step.eq.(iproc-1)) write(Invar%stdlog,*) MPIdata%me_step,data_gather(:,:) 234 !FB call flush_unit(6) 235 !FB call xmpi_barrier(MPIdata%comm_step) 236 !FB end do 237 !FB write(Invar%stdlog,*) '====== END =======' 238 !FB call flush_unit(6) 239 !FB call abi_abort("COLL") 240 !TEST 241 242 if (args%dry_run /= 0) then 243 call wrtout(std_out, "Dry run mode. Exiting after have read the input") 244 goto 100 245 end if 246 247 ! Initialize basic quantities 248 print_mem_report = 1 249 natom = Invar%natom 250 natom_unitcell = Invar%natom_unitcell 251 stdout = Invar%stdout 252 stdlog = Invar%stdlog 253 nshell_max = 500 254 255 !========================================================================================== 256 !============== Define the ideal lattice, symmetries and Brillouin zone =================== 257 !========================================================================================== 258 ! Define all the quantities needed to buid the lattice (rprim*, acell*, brav*...) 259 call tdep_make_latt(Invar,Lattice) 260 261 ! Compute all the symmetries coming from the bravais lattice 262 call tdep_make_sym(Invar,Lattice,MPIdata,Sym) 263 264 ! Initialize the Brillouin zone and compute the q-points path 265 call tdep_make_qptpath(Invar,Lattice,MPIdata,Qpt) 266 267 !========================================================================================== 268 !======== 1/ Determine ideal positions and distances ====================================== 269 !======== 2/ Find the matching between the ideal and average ============================== 270 !======== (from the MD simulations) positions. ========================================== 271 !======== 3/ Find the symmetry operation between the reference and image bonds ============ 272 !======== 4/ Write output quantities needed to visualize the neighbouring distances ======= 273 !========================================================================================== 274 ABI_MALLOC(Rlatt4Abi ,(3,natom_unitcell,natom)) ; Rlatt4Abi (:,:,:)=0.d0 275 ABI_MALLOC(distance,(natom,natom,4)) ; distance(:,:,:)=0.d0 276 ABI_MALLOC(Rlatt_cart,(3,natom_unitcell,natom)) ; Rlatt_cart(:,:,:)=0.d0 277 ABI_MALLOC(ucart,(3,natom,Invar%my_nstep)) ; ucart(:,:,:)=0.d0 278 ABI_MALLOC(Forces_MD,(3*natom*Invar%my_nstep)) ; Forces_MD(:)=0.d0 279 280 call tdep_MatchIdeal2Average(distance,Forces_MD,Invar,Lattice,MPIdata,Rlatt_cart,Rlatt4Abi,Sym,ucart) 281 282 !========================================================================================== 283 !============== Initialize Crystal and DDB ABINIT Datatypes =============================== 284 !========================================================================================== 285 call tdep_init_crystal(Crystal,Invar,Lattice,Sym) 286 call tdep_init_ddb(Crystal,DDB,Invar,Lattice,MPIdata,Qbz) 287 288 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 289 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 290 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=# CALCULATION OF THE 2nd ORDER =#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 291 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 292 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 293 294 !========================================================================================== 295 !============== Initialize the Shell1at datatype ========================================== 296 !========================================================================================== 297 ABI_MALLOC(proj_tmp,(3,3,nshell_max)) ; proj_tmp(:,:,:)=0.d0 298 call tdep_init_shell1at(distance,Invar,MPIdata,3,nshell_max,ncoeff1st,1,proj_tmp,Shell1at,Sym) 299 ABI_MALLOC(proj1st ,(3,3,Shell1at%nshell)) ; proj1st(:,:,:)=0.d0 300 do ishell=1,Shell1at%nshell 301 proj1st(:,:,ishell) = proj_tmp(:,:,ishell) 302 end do 303 ABI_FREE(proj_tmp) 304 !Rotational invariances (1st order) 305 ! constraints = 3 306 CoeffMoore%nconst_1st=3**2 307 308 !========================================================================================== 309 !============== Initialize the Shell2at datatype ========================================== 310 !========================================================================================== 311 ABI_MALLOC(proj_tmp,(9,9,nshell_max)) ; proj_tmp(:,:,:)=0.d0 312 call tdep_init_shell2at(distance,Invar,MPIdata,9,nshell_max,ncoeff2nd,2,proj_tmp,Shell2at,Sym) 313 ABI_MALLOC(proj2nd ,(9,9,Shell2at%nshell)) ; proj2nd(:,:,:)=0.d0 314 do ishell=1,Shell2at%nshell 315 proj2nd(:,:,ishell) = proj_tmp(:,:,ishell) 316 end do 317 ABI_FREE(proj_tmp) 318 !Rotational invariances (2nd order) + Symetry of the Dynamical Matrix + Huang invariances 319 ! constraints = natom*3**2 + (3*natom_unitcell)**2 + 3**4 320 CoeffMoore%nconst_rot2nd = 3**3*natom_unitcell 321 CoeffMoore%nconst_dynmat = (3*natom_unitcell)**2 322 CoeffMoore%nconst_huang = 3**4 323 CoeffMoore%nconst_2nd = CoeffMoore%nconst_rot2nd + CoeffMoore%nconst_dynmat + CoeffMoore%nconst_huang 324 325 !========================================================================================== 326 !============== Initialize the IFC Abinit datatype ======================================== 327 !========================================================================================== 328 ABI_MALLOC(Phi1,(3*natom)); Phi1(:) =0.d0 329 call tdep_init_phi2(Phi2,Invar%loto,natom) 330 call tdep_init_ifc(Crystal,DDB,Ifc,Invar,Lattice,MPIdata,Phi2,Rlatt4Abi,Shell2at,Sym) 331 332 !========================================================================================== 333 !============== Initialize the Shell3at datatype ========================================== 334 !========================================================================================== 335 ncoeff3rd=0 336 CoeffMoore%nconst_3rd=0 337 Shell3at%nshell=1 338 if (Invar%order.ge.3) then 339 ABI_MALLOC(proj_tmp,(27,27,nshell_max)) ; proj_tmp(:,:,:)=0.d0 340 call tdep_init_shell3at(distance,Invar,MPIdata,27,nshell_max,ncoeff3rd,3,proj_tmp,Shell3at,Sym) 341 ABI_MALLOC(proj3rd ,(27,27,Shell3at%nshell)) ; proj3rd(:,:,:)=0.d0 342 do ishell=1,Shell3at%nshell 343 proj3rd(:,:,ishell) = proj_tmp(:,:,ishell) 344 end do 345 ABI_FREE(proj_tmp) 346 ! Rotational invariances (3rd order) + acoustic sum rules (3rd order) 347 ! constraints = natom_unitcell*natom*3**4 + natom_unitcell*natom*3**3 348 CoeffMoore%nconst_rot3rd = 3**4*natom_unitcell*natom 349 CoeffMoore%nconst_asr3rd = 3**3*natom_unitcell*natom 350 CoeffMoore%nconst_3rd = CoeffMoore%nconst_rot3rd + CoeffMoore%nconst_asr3rd 351 end if 352 353 !========================================================================================== 354 !============== Initialize the Shell4at datatype ========================================== 355 !========================================================================================== 356 ncoeff4th=0 357 CoeffMoore%nconst_4th=0 358 Shell4at%nshell=1 359 if (Invar%order==4) then 360 ABI_MALLOC(proj_tmp,(81,81,nshell_max)) ; proj_tmp(:,:,:)=0.d0 361 call tdep_init_shell4at(distance,Invar,MPIdata,81,nshell_max,ncoeff4th,4,proj_tmp,Shell4at,Sym) 362 ABI_MALLOC(proj4th ,(81,81,Shell4at%nshell)) ; proj4th(:,:,:)=0.d0 363 do ishell=1,Shell4at%nshell 364 proj4th(:,:,ishell) = proj_tmp(:,:,ishell) 365 end do 366 ABI_FREE(proj_tmp) 367 ! Rotational invariances (4th order) + acoustic sum rules (4th order) 368 ! constraints = natom_unitcell*natom**2*3**5 + natom_unitcell*natom**2*3**4 369 !FB CoeffMoore%nconst_rot4th = 3**5*natom_unitcell*natom**2 370 !FB4TH CoeffMoore%nconst_asr4th = 3**4*natom_unitcell*natom**2 371 CoeffMoore%nconst_rot4th = 0 372 CoeffMoore%nconst_asr4th = 0 373 CoeffMoore%nconst_4th = CoeffMoore%nconst_rot4th + CoeffMoore%nconst_asr4th 374 end if 375 CoeffMoore%ncoeff1st=ncoeff1st 376 CoeffMoore%ncoeff2nd=ncoeff2nd 377 CoeffMoore%ncoeff3rd=ncoeff3rd 378 CoeffMoore%ncoeff4th=ncoeff4th 379 380 !========================================================================================== 381 !================= Build fcoeff and compute constraints =================================== 382 !========================================================================================== 383 ABI_MALLOC(Forces_TDEP,(3*Invar%natom*Invar%my_nstep)); Forces_TDEP(:)=0.d0 384 ABI_MALLOC(Fresid ,(3*Invar%natom*Invar%my_nstep)); Fresid(:)=Forces_MD(:) 385 ABI_MALLOC(Phi1Ui ,(Invar%my_nstep)); Phi1Ui (:)=0.d0 386 ABI_MALLOC(Phi2UiUj ,(Invar%my_nstep)); Phi2UiUj (:)=0.d0 387 ABI_MALLOC(Phi3UiUjUk ,(Invar%my_nstep)); Phi3UiUjUk (:)=0.d0 388 ABI_MALLOC(Phi4UiUjUkUl,(Invar%my_nstep)); Phi4UiUjUkUl(:)=0.d0 389 ntotcoeff=CoeffMoore%ncoeff1st + CoeffMoore%ncoeff2nd + CoeffMoore%ncoeff3rd + CoeffMoore%ncoeff4th 390 ntotconst=CoeffMoore%nconst_1st + CoeffMoore%nconst_2nd + CoeffMoore%nconst_3rd + CoeffMoore%nconst_4th 391 CoeffMoore%ntotcoeff=ntotcoeff 392 CoeffMoore%ntotconst=ntotconst 393 394 !LOTO 395 !Remove the supercell contribution (the "LR part") included in total forces 396 !before computing the "SR part". The full "LR part" will be added later. 397 if (Invar%loto) then 398 call tdep_ifc2phi2(Ifc%dipdip,Ifc,Invar,Lattice,Invar%natom_unitcell,1,Phi2,Rlatt_cart,Shell2at,Sym) 399 call tdep_calc_ftot2(Forces_TDEP,Invar,Phi1,Phi1Ui,Phi2%Tot,Phi2UiUj,ucart) 400 Fresid(:)=Fresid(:)-Forces_TDEP(:) 401 Phi1 (:)=0.d0 402 Phi1Ui (:)=0.d0 403 Phi2%SR (:,:)=0.d0 404 if (Invar%loto) then 405 Phi2%LR (:,:)=0.d0 406 Phi2%Tot (:,:)=0.d0 407 end if 408 Phi2UiUj (:)=0.d0 409 Forces_TDEP(:)=0.d0 410 end if 411 !LOTO 412 413 do istep=1,Invar%my_nstep 414 do iatom=1,Invar%natom 415 do ii=1,3 416 Fresid(ii+3*(iatom-1)+3*Invar%natom*(istep-1))=Fresid(ii+3*(iatom-1)+3*Invar%natom*(istep-1))*& 417 & Invar%weights(istep) 418 end do 419 end do 420 end do 421 422 ABI_MALLOC(CoeffMoore%fcoeff,(3*natom*Invar%my_nstep,ntotcoeff)); CoeffMoore%fcoeff(:,:)=0.d0 423 if (Invar%readifc.ne.1) then 424 call tdep_calc_phi1fcoeff(CoeffMoore,Invar,proj1st,Shell1at,Sym) 425 call tdep_calc_phi2fcoeff(CoeffMoore,Invar,proj2nd,Shell2at,Sym,ucart) 426 end if 427 if (Invar%order.ge.3) then 428 ABI_MALLOC(Phi3_ref,(3,3,3,Shell3at%nshell)) ; Phi3_ref(:,:,:,:)=0.d0 429 call tdep_calc_phi3fcoeff(CoeffMoore,Invar,proj3rd,Shell3at,Sym,ucart) 430 end if 431 if (Invar%order.eq.4) then 432 ABI_MALLOC(Phi4_ref,(3,3,3,3,Shell4at%nshell)); Phi4_ref(:,:,:,:,:)=0.d0 433 call tdep_calc_phi4fcoeff(CoeffMoore,Invar,proj4th,Shell4at,Sym,ucart) 434 end if 435 if (Invar%order.eq.2) then 436 call tdep_calc_constraints(CoeffMoore,distance,Invar,MPIdata,Shell1at%nshell,Shell2at%nshell,& 437 & Shell3at%nshell,Shell4at%nshell,proj1st,proj2nd,& 438 & Shell1at,Shell2at,Sym) 439 else if (Invar%order.eq.3) then 440 call tdep_calc_constraints(CoeffMoore,distance,Invar,MPIdata,Shell1at%nshell,Shell2at%nshell,& 441 & Shell3at%nshell,Shell4at%nshell,proj1st,proj2nd,& 442 & Shell1at,Shell2at,Sym,proj3rd,Shell3at) 443 else if (Invar%order.eq.4) then 444 call tdep_calc_constraints(CoeffMoore,distance,Invar,MPIdata,Shell1at%nshell,Shell2at%nshell,& 445 & Shell3at%nshell,Shell4at%nshell,proj1st,proj2nd,& 446 & Shell1at,Shell2at,Sym,proj3rd,Shell3at,proj4th,Shell4at) 447 end if 448 449 !========================================================================================== 450 !============= Compute the pseudo inverse using the Moore-Penrose method ================== 451 !========================================================================================== 452 ABI_MALLOC(MP_coeff,(ntotcoeff,1)); MP_coeff(:,:)=0.d0 453 !========================================================================================== 454 !=================== If all the Orders are solved simultaneously ========================== 455 !========================================================================================== 456 if (Invar%together.eq.1) then 457 write(stdout,*) '############### (Solve simultaneously all the orders) #######################' 458 if (Invar%readifc.ne.1) then 459 call tdep_calc_MoorePenrose(CoeffMoore,Fresid,0,Invar,MP_coeff,MPIdata) 460 end if 461 ABI_MALLOC(Phi1_coeff,(ncoeff1st,1)); Phi1_coeff(:,:)=MP_coeff(1:ncoeff1st,:) 462 ABI_MALLOC(Phi2_coeff,(ncoeff2nd,1)); Phi2_coeff(:,:)=MP_coeff(ncoeff1st+1:ncoeff1st+ncoeff2nd,:) 463 if (Invar%readifc.ne.1) then 464 call tdep_calc_phi1(Invar,ncoeff1st,proj1st,Phi1_coeff,Phi1 ,Shell1at,Sym) 465 call tdep_calc_phi2(Invar,ncoeff2nd,proj2nd,Phi2_coeff,Phi2%SR,Shell2at,Sym) 466 end if 467 ABI_FREE(proj1st) 468 ABI_FREE(proj2nd) 469 ABI_FREE(Phi1_coeff) 470 ABI_FREE(Phi2_coeff) 471 call tdep_calc_ftot2(Forces_TDEP,Invar,Phi1,Phi1Ui,Phi2%SR,Phi2UiUj,ucart) 472 if (Invar%order.ge.3) then 473 ABI_MALLOC(Phi3_coeff,(ncoeff3rd,1)); Phi3_coeff(:,:)=0.d0 474 Phi3_coeff(:,:)=MP_coeff(ncoeff1st+ncoeff2nd+1:ncoeff1st+ncoeff2nd+ncoeff3rd,:) 475 call tdep_calc_phi3ref(ncoeff3rd,proj3rd,Phi3_coeff,Phi3_ref,Shell3at) 476 ABI_FREE(proj3rd) 477 ABI_FREE(Phi3_coeff) 478 call tdep_calc_ftot3(Forces_TDEP,Invar,Phi3_ref,Phi3UiUjUk,Shell3at,ucart,Sym) 479 end if 480 if (Invar%order.ge.4) then 481 ABI_MALLOC(Phi4_coeff,(ncoeff4th,1)); Phi4_coeff(:,:)=0.d0 482 Phi4_coeff(:,:)=MP_coeff(ncoeff1st+ncoeff2nd+ncoeff3rd+1:ntotcoeff,:) 483 call tdep_calc_phi4ref(ncoeff4th,proj4th,Phi4_coeff,Phi4_ref,Shell4at) 484 ABI_FREE(proj4th) 485 ABI_FREE(Phi4_coeff) 486 call tdep_calc_ftot4(Forces_TDEP,Invar,Phi4_ref,Phi4UiUjUkUl,Shell4at,ucart,Sym) 487 end if 488 489 !========================================================================================== 490 !=================== If all the Orders are solved successively ============================ 491 !========================================================================================== 492 !ATTENTION : Le LR semble enleve a l'ordre 2 mais pas a l'ordre 3 et 4. On repart de 493 ! Forces_MD tout en bas et pas de Fresid (car les ordres 2 et 3 sont supprimes au 494 else if (Invar%together.eq.0) then 495 write(stdout,*) '################## (Solve successively each order) ##########################' 496 do ii=1,Invar%order-1 497 write(stdout,*) ' For order=',ii+1 498 if (Invar%readifc.eq.1) cycle 499 call tdep_calc_MoorePenrose(CoeffMoore,Fresid,ii,Invar,MP_coeff,MPIdata) 500 if (ii.eq.1) then 501 ABI_MALLOC(Phi1_coeff,(ncoeff1st,1)); Phi1_coeff(:,:)=MP_coeff(1:ncoeff1st,:) 502 ABI_MALLOC(Phi2_coeff,(ncoeff2nd,1)); Phi2_coeff(:,:)=MP_coeff(ncoeff1st+1:ncoeff1st+ncoeff2nd,:) 503 if (Invar%readifc.ne.1) then 504 call tdep_calc_phi1(Invar,ncoeff1st,proj1st,Phi1_coeff,Phi1 ,Shell1at,Sym) 505 call tdep_calc_phi2(Invar,ncoeff2nd,proj2nd,Phi2_coeff,Phi2%SR,Shell2at,Sym) 506 end if 507 ABI_FREE(proj1st) 508 ABI_FREE(proj2nd) 509 ABI_FREE(Phi1_coeff) 510 ABI_FREE(Phi2_coeff) 511 call tdep_calc_ftot2(Forces_TDEP,Invar,Phi1,Phi1Ui,Phi2%SR,Phi2UiUj,ucart) 512 else if (ii.eq.2) then 513 ABI_MALLOC(Phi3_coeff,(ncoeff3rd,1)); Phi3_coeff(:,:)=0.d0 514 Phi3_coeff(:,:)=MP_coeff(ncoeff1st+ncoeff2nd+1:ncoeff1st+ncoeff2nd+ncoeff3rd,:) 515 call tdep_calc_phi3ref(ncoeff3rd,proj3rd,Phi3_coeff,Phi3_ref,Shell3at) 516 ABI_FREE(proj3rd) 517 ABI_FREE(Phi3_coeff) 518 call tdep_calc_ftot3(Forces_TDEP,Invar,Phi3_ref,Phi3UiUjUk,Shell3at,ucart,Sym) 519 else if (ii.eq.3) then 520 ABI_MALLOC(Phi4_coeff,(ncoeff4th,1)); Phi4_coeff(:,:)=0.d0 521 Phi4_coeff(:,:)=MP_coeff(ncoeff1st+ncoeff2nd+ncoeff3rd+1:ntotcoeff,:) 522 call tdep_calc_phi4ref(ncoeff4th,proj4th,Phi4_coeff,Phi4_ref,Shell4at) 523 ABI_FREE(proj4th) 524 ABI_FREE(Phi4_coeff) 525 call tdep_calc_ftot4(Forces_TDEP,Invar,Phi4_ref,Phi4UiUjUkUl,Shell4at,ucart,Sym) 526 end if 527 do istep=1,Invar%my_nstep 528 do iatom=1,Invar%natom 529 do jj=1,3 530 Fresid(jj+3*(iatom-1)+3*Invar%natom*(istep-1))=(Forces_MD(jj+3*(iatom-1)+3*Invar%natom*(istep-1)) -& 531 & Forces_TDEP(jj+3*(iatom-1)+3*Invar%natom*(istep-1)))*& 532 & Invar%weights(istep) 533 end do ! ii 534 end do ! iatom 535 end do ! istep 536 end do 537 end if 538 ABI_FREE(CoeffMoore%const) 539 ABI_FREE(CoeffMoore%fcoeff) 540 ABI_FREE(Fresid) 541 ABI_FREE(MP_coeff) 542 ABI_FREE(ucart) 543 call tdep_destroy_shell(natom,1,Shell1at) 544 545 !LOTO 546 if (Invar%loto) then 547 Phi2%Tot=Phi2%LR+Phi2%SR 548 end if 549 !LOTO 550 551 !========================================================================================== 552 !=================== Write the IFC and check the constraints ============================== 553 !========================================================================================== 554 call tdep_write_phi1(Invar,Phi1) 555 call tdep_write_phi2(distance,Invar,MPIdata,Phi2%SR,Shell2at) 556 if (Invar%order.ge.3) then 557 call tdep_write_phi3(distance,Invar,Phi3_ref,Shell3at,Sym) 558 end if 559 if (Invar%order.ge.4) then 560 call tdep_write_phi4(distance,Invar,Phi4_ref,Shell4at,Sym) 561 end if 562 563 if (Invar%order.eq.2) then 564 call tdep_check_constraints(distance,Invar,Phi2%SR,Phi1,Shell3at%nshell,Shell4at%nshell,Sym) 565 else if (Invar%order.eq.3) then 566 call tdep_check_constraints(distance,Invar,Phi2%SR,Phi1,Shell3at%nshell,Shell4at%nshell,Sym,& 567 & Phi3_ref,Shell3at) 568 else if (Invar%order.eq.4) then 569 call tdep_check_constraints(distance,Invar,Phi2%SR,Phi1,Shell3at%nshell,Shell4at%nshell,Sym,& 570 & Phi3_ref,Shell3at,Phi4_ref,Shell4at) 571 end if 572 ABI_FREE(Phi1) 573 574 !========================================================================================== 575 !===================== Compute the phonon spectrum, the DOS, ============================== 576 !===================== the dynamical matrix and write them =============================== 577 !========================================================================================== 578 write(stdout,*)' ' 579 write(stdout,*) '#############################################################################' 580 write(stdout,*) '############## Compute the phonon spectrum, the DOS, ########################' 581 write(stdout,*) '############## the dynamical matrix and write them ########################' 582 write(stdout,*) '#############################################################################' 583 call tdep_init_eigen2nd(Eigen2nd_MP,Invar%natom_unitcell,Qbz%nqbz) 584 !FB call tdep_init_eigen2nd(Eigen2nd_MP,Invar%natom_unitcell,Qbz%nqibz) 585 call tdep_init_eigen2nd(Eigen2nd_path,Invar%natom_unitcell,Qpt%nqpt) 586 call tdep_calc_phdos(Crystal,DDB,Eigen2nd_MP,Eigen2nd_path,Ifc,Invar,Lattice,MPIdata,natom,& 587 & natom_unitcell,Phi2,PHdos,Qbz,Qpt,Rlatt4abi,Shell2at,Sym) 588 call tdep_destroy_shell(natom,2,Shell2at) 589 ABI_FREE(Rlatt4Abi) 590 write(stdout,'(a)') ' See the dij.dat, omega.dat and eigenvectors files' 591 write(stdout,'(a)') ' See also the DDB file' 592 593 !========================================================================================== 594 !===================== Compute the elastic constants ====================================== 595 !========================================================================================== 596 call tdep_calc_elastic(Phi2%SR,distance,Invar,Lattice) 597 call tdep_destroy_phi2(Phi2,Invar%loto) 598 599 !========================================================================================== 600 !=========== Compute U_0, the "free energy" and the forces (from the model) =============== 601 !========================================================================================== 602 call tdep_calc_model(Forces_MD,Forces_TDEP,Invar,MPIdata,Phi1Ui,Phi2UiUj,& 603 & Phi3UiUjUk,Phi4UiUjUkUl,U0) 604 ABI_FREE(Forces_MD) 605 ABI_FREE(Forces_TDEP) 606 ABI_FREE(Phi1Ui) 607 ABI_FREE(Phi2UiUj) 608 ABI_FREE(Phi3UiUjUk) 609 ABI_FREE(Phi4UiUjUkUl) 610 611 !========================================================================================== 612 !===================== Compute the thermodynamical quantities ============================= 613 !========================================================================================== 614 call tdep_calc_thermo(Invar,Lattice,MPIdata,PHdos,U0) 615 call PHdos%free() 616 617 618 if (Invar%order==2) then 619 620 ABI_FREE(distance) 621 ABI_FREE(Rlatt_cart) 622 call Ifc%free() 623 call crystal_free(Crystal) 624 call tdep_destroy_eigen2nd(Eigen2nd_path) 625 call tdep_destroy_eigen2nd(Eigen2nd_MP) 626 call tdep_destroy_sym(Sym) 627 call tdep_destroy_qbz(Qbz) 628 call tdep_destroy_qpt(Qpt) 629 call tdep_destroy_ddb(DDB) 630 call tdep_destroy_invar(Invar) 631 call tdep_destroy_mpidata(MPIdata) 632 633 call tdep_print_Aknowledgments(Invar) 634 call delete_file(Invar%foo,ierr) 635 call delete_file('fort.8',ierr) 636 call flush_unit(stdout) 637 close(unit=stdout) 638 call abinit_doctor(trim(Invar%output_prefix), print_mem_report=print_mem_report) 639 call flush_unit(stdlog) 640 call xmpi_end() 641 stop 642 end if 643 644 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 645 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 646 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=# CALCULATION OF THE 3rd ORDER =#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 647 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 648 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 649 650 if (MPIdata%iam_master) then 651 call tdep_write_gruneisen(distance,Eigen2nd_path,Invar,Phi3_ref,Qpt,Rlatt_cart,Shell3at,Sym) 652 end if 653 call tdep_calc_alpha_gamma(distance,Eigen2nd_MP,Invar,Lattice,MPIdata,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym) 654 655 !FB Begin Lifetime 656 !FB call tdep_calc_lifetime1(Crystal,distance,Eigen2nd_MP,Ifc,Invar,Lattice,Phi3_ref,Qbz,Rlatt_cart,Shell3at,Sym) 657 !FB End Lifetime 658 659 ABI_FREE(Rlatt_cart) 660 ABI_FREE(distance) 661 call tdep_destroy_eigen2nd(Eigen2nd_path) 662 call tdep_destroy_eigen2nd(Eigen2nd_MP) 663 call tdep_destroy_shell(natom,3,Shell3at) 664 ABI_FREE(Phi3_ref) 665 if (Invar%order.eq.4) then 666 call tdep_destroy_shell(natom,4,Shell4at) 667 ABI_FREE(Phi4_ref) 668 end if 669 call Ifc%free() 670 call crystal_free(Crystal) 671 call tdep_destroy_sym(Sym) 672 call tdep_destroy_qbz(Qbz) 673 call tdep_destroy_qpt(Qpt) 674 call tdep_destroy_ddb(DDB) 675 call tdep_destroy_invar(Invar) 676 call tdep_destroy_mpidata(MPIdata) 677 678 !========================================================================================== 679 !================= Write the last informations (aknowledgments...) ======================= 680 !========================================================================================== 681 call tdep_print_Aknowledgments(Invar) 682 call delete_file(Invar%foo,ierr) 683 call delete_file('fort.8',ierr) 684 call flush_unit(stdout) 685 close(unit=stdout) 686 call abinit_doctor(trim(Invar%output_prefix), print_mem_report=print_mem_report) 687 call flush_unit(stdlog) 688 100 call xmpi_end() 689 690 end program atdep