TABLE OF CONTENTS
ABINIT/tdep [ Programs ]
NAME
tdep
FUNCTION
Calculations of phonons using molecular dynamic simulations
COPYRIGHT
Copyright (C) 1998-2020 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)
PARENTS
CHILDREN
tdep_build_phijnn,tdep_calc_psijtot,tdep_calc_alpha_gamma
tdep_calc_constraints,tdep_calc_dij,tdep_calc_elastic,tdep_calc_model
tdep_calc_moorepenrose,tdep_calc_phdos,tdep_calc_phijfcoeff
tdep_calc_psijfcoeff,tdep_calc_thermo,tdep_destroy_shell
tdep_destroy_eigen2nd,tdep_init_crystal,tdep_init_ddb
tdep_init_eigen2nd,tdep_init_ifc,tdep_init_shell2at
tdep_init_shell3at,tdep_make_latt,tdep_make_qptpath
tdep_make_sym,tdep_matchideal2average,tdep_print_aknowledgments
tdep_readecho,tdep_write_gruneisen,tdep_write_dij,tdep_write_yaml
xmpi_end,xmpi_init
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 program tdep 52 53 use defs_basis 54 use m_abicore 55 use m_phonons 56 use m_xmpi 57 use m_io_tools 58 use m_errors 59 use m_argparse 60 61 use m_ifc, only : ifc_type 62 use m_crystal, only : crystal_t 63 use m_ddb, only : ddb_type 64 use m_dynmat, only : ftifc_r2q, ftifc_q2r, asrif9 65 use m_tdep_abitypes, only : tdep_init_crystal, tdep_init_ifc, tdep_init_ddb, tdep_write_ifc 66 use m_tdep_psij, only : tdep_calc_psijfcoeff, tdep_calc_psijtot, tdep_calc_alpha_gamma, tdep_write_gruneisen 67 use m_tdep_phij, only : tdep_calc_phijfcoeff, tdep_calc_pijfcoeff, tdep_build_phijNN, tdep_calc_dij, tdep_write_dij, & 68 Eigen_Variables_type, tdep_init_eigen2nd, tdep_destroy_eigen2nd, & 69 tdep_write_yaml, tdep_build_pijN 70 use m_tdep_latt, only : tdep_make_latt, Lattice_Variables_type 71 use m_tdep_sym, only : tdep_make_sym, Symetries_Variables_type 72 use m_tdep_readwrite, only : tdep_print_Aknowledgments, tdep_ReadEcho, Input_Variables_type 73 use m_tdep_utils, only : Coeff_Moore_type, tdep_calc_MoorePenrose, tdep_MatchIdeal2Average, tdep_calc_model 74 use m_tdep_qpt, only : tdep_make_qptpath, Qpoints_type 75 use m_tdep_phdos, only : tdep_calc_phdos,tdep_calc_elastic,tdep_calc_thermo 76 use m_tdep_shell, only : Shell_Variables_type, tdep_init_shell2at, tdep_init_shell3at, tdep_init_shell1at, & 77 tdep_destroy_shell 78 use m_tdep_constraints, only : tdep_calc_constraints, tdep_check_constraints 79 80 implicit none 81 82 integer :: natom,natom_unitcell,ncoeff1st,ncoeff2nd,ncoeff3rd,ntotcoeff,ntotconst 83 integer :: stdout,iqpt,nshell_max 84 double precision :: U0,Free_Anh 85 double precision, allocatable :: ucart(:,:,:),proj1st(:,:,:),proj2nd(:,:,:),proj3rd(:,:,:) 86 double precision, allocatable :: proj_tmp(:,:,:),Forces_TDEP(:) 87 !FB double precision, allocatable :: fcoeff(:,:),Phij_coeff(:,:),Forces_MD(:),Phij_NN(:,:) 88 double precision, allocatable :: Phij_coeff(:,:),Forces_MD(:),Phij_NN(:,:),Pij_N(:),Pij_coeff(:,:) 89 double precision, allocatable :: Psij_coeff(:,:),Psij_ref(:,:,:,:),MP_coeff(:,:) 90 double precision, allocatable :: distance(:,:,:),Rlatt_cart(:,:,:),Rlatt4Abi(:,:,:) 91 double precision, allocatable :: omega (:),ftot3(:,:) 92 double precision :: qpt_cart(3) 93 double complex , allocatable :: dij(:,:),eigenV(:,:) 94 type(args_t) :: args 95 type(phonon_dos_type) :: PHdos 96 type(Input_Variables_type) :: InVar 97 type(Lattice_Variables_type) :: Lattice 98 type(Symetries_Variables_type) :: Sym 99 type(Qpoints_type) :: Qpt 100 type(ifc_type) :: Ifc 101 type(ddb_type) :: DDB 102 type(crystal_t) :: Crystal 103 type(Shell_Variables_type) :: Shell1at,Shell2at,Shell3at 104 type(Coeff_Moore_type) :: CoeffMoore 105 type(Eigen_Variables_type) :: Eigen2nd 106 107 !****************************************************************** 108 109 !========================================================================================== 110 !===================== Initialization & Reading ========================================== 111 !========================================================================================== 112 call xmpi_init() 113 114 ! Parse command line arguments. 115 args = args_parser(); if (args%exit /= 0) goto 100 116 117 ! Initialize memory profiling if activated at configure time. 118 ! if a full report is desired, set the argument of abimem_init to "2" instead of "0" via the command line. 119 ! note that the file can easily be multiple GB in size so don't use this option normally 120 #ifdef HAVE_MEM_PROFILING 121 call abimem_init(args%abimem_level, limit_mb=args%abimem_limit_mb) 122 #endif 123 124 ! Read input values from the input.in input file 125 call tdep_ReadEcho(InVar) 126 127 if (args%dry_run /= 0) then 128 call wrtout(std_out, "Dry run mode. Exiting after have read the input") 129 goto 100 130 end if 131 132 ! Initialize basic quantities 133 natom =InVar%natom 134 natom_unitcell=InVar%natom_unitcell 135 stdout =InVar%stdout 136 nshell_max =500 137 138 !========================================================================================== 139 !============== Define the ideal lattice, symmetries and Brillouin zone =================== 140 !========================================================================================== 141 ! Define all the quantities needed to buid the lattice (rprim*, acell*, brav*...) 142 call tdep_make_latt(InVar,Lattice) 143 144 ! Compute all the symmetries coming from the bravais lattice 145 call tdep_make_sym(Invar,Lattice,Sym) 146 147 ! Initialize the Brillouin zone and compute the q-points path 148 call tdep_make_qptpath(InVar,Lattice,Qpt) 149 150 !========================================================================================== 151 !======== 1/ Determine ideal positions and distances ====================================== 152 !======== 2/ Find the matching between the ideal and average ============================== 153 !======== (from the MD simulations) positions. ========================================== 154 !======== 3/ Find the symmetry operation between the reference and image bonds ============= 155 !======== 4/ Write output quantities needed to visualize the neighbouring distances ======= 156 !========================================================================================== 157 ABI_MALLOC(Rlatt4Abi ,(3,natom_unitcell,natom)) ; Rlatt4Abi (:,:,:)=0.d0 158 ABI_MALLOC(distance,(natom,natom,4)) ; distance(:,:,:)=0.d0 159 ABI_MALLOC(Rlatt_cart,(3,natom_unitcell,natom)) ; Rlatt_cart(:,:,:)=0.d0 160 ABI_MALLOC(ucart,(3,natom,InVar%nstep)) ; ucart(:,:,:)=0.d0 161 ABI_MALLOC(Forces_MD,(3*natom*InVar%nstep)) ; Forces_MD(:)=0.d0 162 163 call tdep_MatchIdeal2Average(distance,Forces_MD,InVar,Lattice,Rlatt_cart,Rlatt4Abi,Sym,ucart) 164 165 !========================================================================================== 166 !============== Initialize Crystal and DDB ABINIT Datatypes =============================== 167 !========================================================================================== 168 call tdep_init_crystal(Crystal,InVar,Lattice,Sym) 169 call tdep_init_ddb(DDB,InVar,Lattice) 170 171 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 172 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 173 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=# CALCULATION OF THE 2nd ORDER =#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 174 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 175 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 176 177 !========================================================================================== 178 !============== Initialize the Shell1at datatype ========================================== 179 !========================================================================================== 180 ABI_MALLOC(proj_tmp,(3,3,nshell_max)) ; proj_tmp(:,:,:)=0.d0 181 call tdep_init_shell1at(distance,InVar,3,nshell_max,ncoeff1st,1,proj_tmp,Shell1at,Sym) 182 ABI_MALLOC(proj1st ,(3,3,Shell1at%nshell)) ; proj1st(:,:,:)=0.d0 183 proj1st = reshape (proj_tmp, (/ 3,3,Shell1at%nshell /)) 184 ABI_FREE(proj_tmp) 185 !Rotational invariances (1st order) 186 ! constraints = 3 187 CoeffMoore%nconst_1st=3 188 ntotconst=CoeffMoore%nconst_1st 189 190 !========================================================================================== 191 !============== Initialize the Shell2at datatype ========================================== 192 !========================================================================================== 193 ABI_MALLOC(proj_tmp,(9,9,nshell_max)) ; proj_tmp(:,:,:)=0.d0 194 call tdep_init_shell2at(distance,InVar,9,nshell_max,ncoeff2nd,2,proj_tmp,Shell2at,Sym) 195 ABI_MALLOC(proj2nd ,(9,9,Shell2at%nshell)) ; proj2nd(:,:,:)=0.d0 196 proj2nd = reshape (proj_tmp, (/ 9,9,Shell2at%nshell /)) 197 ABI_FREE(proj_tmp) 198 !Rotational invariances (2nd order) + Symetry of the Dynamical Matrix + Huang invariances 199 ! constraints = natom*3**2 + (3*natom_unitcell)**2 + 3**4 200 CoeffMoore%nconst_rot2nd = natom_unitcell*9 201 CoeffMoore%nconst_dynmat = (3*natom_unitcell)**2 202 CoeffMoore%nconst_huang = 81 203 CoeffMoore%nconst_2nd = CoeffMoore%nconst_rot2nd + CoeffMoore%nconst_dynmat + CoeffMoore%nconst_huang 204 ntotconst=CoeffMoore%nconst_1st + CoeffMoore%nconst_2nd 205 206 !========================================================================================== 207 !============== Initialize the IFC Abinit datatype ======================================== 208 !========================================================================================== 209 ABI_MALLOC(Pij_N ,(3*natom)) ; Pij_N (:) =0.d0 210 ABI_MALLOC(Phij_NN,(3*natom,3*natom)) ; Phij_NN(:,:)=0.d0 211 call tdep_init_ifc(Crystal,DDB,Ifc,InVar,Lattice,Phij_NN,Rlatt4Abi,Shell2at,Sym) 212 213 !========================================================================================== 214 !============== Initialize the Shell3at datatype (if SC_order==1) ========================= 215 !========================================================================================== 216 ncoeff3rd=0 217 if (InVar%Order==3) then 218 ABI_MALLOC(proj_tmp,(27,27,nshell_max)) ; proj_tmp(:,:,:)=0.d0 219 call tdep_init_shell3at(distance,InVar,27,nshell_max,ncoeff3rd,3,proj_tmp,Shell3at,Sym) 220 ABI_MALLOC(proj3rd ,(27,27,Shell3at%nshell)) ; proj3rd(:,:,:)=0.d0 221 proj3rd = reshape (proj_tmp, (/ 27,27,Shell3at%nshell /)) 222 ABI_FREE(proj_tmp) 223 ! Rotational invariances (3rd order) + acoustic sum rules (3rd order) 224 ! constraints = natom_unitcell*natom*3**3) + 3permutations*3**3*nshell2at 225 !FB CoeffMoore%nconst_rot3rd = natom_unitcell*natom*3**3 226 CoeffMoore%nconst_rot3rd = 3*natom_unitcell*natom*3**4 227 CoeffMoore%nconst_asr3rd = 8*natom_unitcell*natom*3**3 228 CoeffMoore%nconst_3rd = CoeffMoore%nconst_rot3rd + CoeffMoore%nconst_asr3rd 229 ntotconst=CoeffMoore%nconst_1st + CoeffMoore%nconst_2nd + CoeffMoore%nconst_3rd 230 end if 231 ntotcoeff=ncoeff1st+ncoeff2nd+ncoeff3rd 232 CoeffMoore%ntotcoeff=ntotcoeff 233 CoeffMoore%ntotconst=ntotconst 234 CoeffMoore%ncoeff1st=ncoeff1st 235 CoeffMoore%ncoeff2nd=ncoeff2nd 236 CoeffMoore%ncoeff3rd=ncoeff3rd 237 ABI_MALLOC(CoeffMoore%fcoeff,(3*natom*InVar%nstep+ntotconst,ntotcoeff)); CoeffMoore%fcoeff(:,:)=0.d0 238 239 !========================================================================================== 240 !============= Build fcoeff, needed for the Moore-Penrose method just below =============== 241 !============================== and compute constraints =================================== 242 !========================================================================================== 243 if (InVar%ReadIFC.ne.1) then 244 call tdep_calc_pijfcoeff(CoeffMoore,InVar,proj1st,Shell1at,Sym) 245 call tdep_calc_phijfcoeff(CoeffMoore,InVar,proj2nd,Shell2at,Sym,ucart) 246 end if 247 if (InVar%Order==3) then 248 call tdep_calc_psijfcoeff(CoeffMoore,InVar,proj3rd,Shell3at,Sym,ucart) 249 CoeffMoore%fcoeff(:,ncoeff2nd+1:ntotcoeff)=CoeffMoore%fcoeff(:,ncoeff2nd+1:ntotcoeff)/2.d0 250 call tdep_calc_constraints(CoeffMoore,distance,InVar,Shell1at%nshell,Shell2at%nshell,Shell3at%nshell,Sym,& 251 & proj1st,Shell1at,proj2nd,Shell2at,proj3rd,Shell3at) 252 else 253 call tdep_calc_constraints(CoeffMoore,distance,InVar,Shell1at%nshell,Shell2at%nshell,Shell3at%nshell,Sym,& 254 & proj1st,Shell1at,proj2nd,Shell2at) 255 end if 256 257 !========================================================================================== 258 !============= Compute the pseudo inverse using the Moore-Penrose method ================== 259 !========================================================================================== 260 ABI_MALLOC(MP_coeff,(ntotcoeff,1)); MP_coeff(:,:)=0.d0 261 if (InVar%ReadIFC.ne.1) then 262 call tdep_calc_MoorePenrose(Forces_MD,CoeffMoore,InVar,MP_coeff) 263 end if 264 ABI_MALLOC(Pij_coeff ,(ncoeff1st,1)); Pij_coeff (:,:)=MP_coeff(1:ncoeff1st,:) 265 ABI_MALLOC(Phij_coeff,(ncoeff2nd,1)); Phij_coeff(:,:)=MP_coeff(ncoeff1st+1:ncoeff1st+ncoeff2nd,:) 266 267 !========================================================================================== 268 !==== Reorganize the IFC coefficients into the whole Pij_N & Phij_NN matrices ============= 269 !=================== and check the constraints ============================================ 270 !========================================================================================== 271 if (InVar%ReadIFC.ne.1) then 272 call tdep_build_pijN(InVar,ncoeff1st,proj1st,Pij_coeff,Pij_N,Shell1at,Sym) 273 call tdep_build_phijNN(distance,InVar,ncoeff2nd,proj2nd,Phij_coeff,Phij_NN,Shell2at,Sym) 274 call tdep_check_constraints(distance,InVar,Phij_NN,Pij_N,1) 275 end if 276 ABI_FREE(proj2nd) 277 ABI_FREE(Pij_coeff) 278 ABI_FREE(Phij_coeff) 279 280 !========================================================================================== 281 !===================== Compute the phonons density of states ============================== 282 !========================================================================================== 283 call tdep_calc_phdos(Crystal,Ifc,InVar,Lattice,natom,natom_unitcell,Phij_NN,PHdos,Qpt,Rlatt4Abi,Shell2at,Sym) 284 call tdep_destroy_shell(natom,2,Shell2at) 285 ABI_FREE(Rlatt4Abi) 286 287 !========================================================================================== 288 !============= Compute the dynamical matrix and phonon spectrum =========================== 289 !================ then print Dij, omega, eigenvectors ===================================== 290 !========================================================================================== 291 write(stdout,*)' ' 292 write(stdout,*) '#############################################################################' 293 write(stdout,*) '######################## Dynamical matrix ###################################' 294 write(stdout,*) '#############################################################################' 295 open(unit=53,file=trim(InVar%output_prefix)//'omega.dat') 296 open(unit=52,file=trim(InVar%output_prefix)//'dij.dat') 297 open(unit=51,file=trim(InVar%output_prefix)//'eigenvectors.dat') 298 ABI_MALLOC(dij ,(3*InVar%natom_unitcell,3*InVar%natom_unitcell)) 299 ABI_MALLOC(eigenV,(3*InVar%natom_unitcell,3*InVar%natom_unitcell)) 300 ABI_MALLOC(omega,(3*InVar%natom_unitcell)) 301 call tdep_init_eigen2nd(Eigen2nd,InVar%natom_unitcell,Qpt%nqpt) 302 do iqpt=1,Qpt%nqpt 303 dij(:,:)=zero ; eigenV(:,:)=zero ; omega(:)=zero 304 qpt_cart(:)=Qpt%qpt_cart(:,iqpt) 305 call tdep_calc_dij (dij,eigenV,iqpt,InVar,Lattice,omega,Phij_NN,qpt_cart,Rlatt_cart) 306 call tdep_write_dij(dij,eigenV,iqpt,InVar,Lattice,omega,qpt_cart) 307 Eigen2nd%eigenval(:,iqpt)= omega(:) 308 Eigen2nd%eigenvec(:,:,iqpt)=eigenV(:,:) 309 end do 310 ABI_FREE(dij) 311 ABI_FREE(eigenV) 312 ABI_FREE(omega) 313 close(53) 314 close(52) 315 close(51) 316 call tdep_write_yaml(Eigen2nd,Qpt,InVar%output_prefix) 317 write(InVar%stdout,'(a)') ' See the dij.dat, omega.dat and eigenvectors files' 318 !========================================================================================== 319 !===================== Compute the elastic constants ====================================== 320 !========================================================================================== 321 call tdep_calc_elastic(Phij_NN,distance,InVar,Lattice) 322 323 !========================================================================================== 324 !=========== Compute U_0, the "free energy" and the forces (from the model) =============== 325 !========================================================================================== 326 ABI_MALLOC(Forces_TDEP,(3*InVar%natom*InVar%nstep)); Forces_TDEP(:)=0.d0 327 call tdep_calc_model(Free_Anh,Forces_MD,Forces_TDEP,InVar,Phij_NN,Pij_N,ucart,U0) 328 329 !========================================================================================== 330 !===================== Compute the thermodynamical quantities ============================= 331 !========================================================================================== 332 call tdep_calc_thermo(Free_Anh,InVar,Lattice,PHdos,U0) 333 334 if (InVar%Order==2) then 335 call tdep_print_Aknowledgments(InVar) 336 call xmpi_end() 337 stop 338 end if 339 340 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 341 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 342 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=# CALCULATION OF THE 3rd ORDER =#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 343 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 344 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#= 345 346 !========================================================================================== 347 !============= Initialize Psij_coeff from Moore-Penrose coefficients ====================== 348 !========================================================================================== 349 ABI_MALLOC(Psij_coeff,(ncoeff3rd,1)); Psij_coeff(:,:)=0.d0 350 Psij_coeff(:,:)=MP_coeff(ncoeff2nd+1:ntotcoeff,:) 351 ABI_FREE(MP_coeff) 352 353 !========================================================================================== 354 !============= Compute all the IFC coefficients of the Psi_ijk matrix ===================== 355 !========================================================================================== 356 ABI_MALLOC(Psij_ref,(3,3,3,Shell3at%nshell)) ; Psij_ref(:,:,:,:)=0.d0 357 call tdep_calc_psijtot(distance,InVar,ncoeff3rd,proj3rd,Psij_coeff,Psij_ref,Shell3at,Sym) 358 ABI_MALLOC(ftot3,(3*InVar%natom,InVar%nstep)); ftot3(:,:)=0.d0 359 call tdep_check_constraints(distance,InVar,Phij_NN,Pij_N,Shell3at%nshell,ftot3,Psij_ref,Shell3at,Sym,ucart) 360 ABI_FREE(Psij_coeff) 361 ABI_FREE(proj3rd) 362 call tdep_write_gruneisen(distance,Eigen2nd,InVar,Lattice,Psij_ref,Qpt,Rlatt_cart,Shell3at,Sym) 363 call tdep_calc_alpha_gamma(Crystal,distance,DDB,Ifc,InVar,Lattice,Psij_ref,Rlatt_cart,Shell3at,Sym) 364 !FB stop 365 ABI_FREE(Rlatt_cart) 366 call tdep_destroy_eigen2nd(Eigen2nd) 367 call tdep_destroy_shell(natom,3,Shell3at) 368 call tdep_calc_model(Free_Anh,Forces_MD,Forces_TDEP,InVar,Phij_NN,Pij_N,ucart,U0,ftot3) 369 ABI_FREE(ftot3) 370 ABI_FREE(distance) 371 ABI_FREE(Forces_MD) 372 ABI_FREE(Forces_TDEP) 373 ABI_FREE(Pij_N) 374 ABI_FREE(Phij_NN) 375 ABI_FREE(Psij_ref) 376 ABI_FREE(ucart) 377 !========================================================================================== 378 !================= Write the last informations (aknowledgments...) ======================= 379 !========================================================================================== 380 call tdep_print_Aknowledgments(InVar) 381 382 call abinit_doctor("__fftprof") 383 384 100 call xmpi_end() 385 386 end program tdep