diff --git a/CMakeLists.txt b/CMakeLists.txt index 45a5aecf..8cb4ea9d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,6 +8,7 @@ cmake_minimum_required(VERSION 3.16) # set cpp standard flags set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_EXTENSIONS OFF) +cmake_policy(SET CMP0144 NEW) project(Phoebe LANGUAGES C CXX Fortran VERSION 1.0) diff --git a/README.md b/README.md index 2aa40a2b..2eeb50fc 100644 --- a/README.md +++ b/README.md @@ -26,21 +26,20 @@ If you feel you've found a bug or seen some unexpected behavior, please let us k ### Current functionalities #### Electronic Transport - * Electron-phonon and phonon-electron scattering by Wannier interpolation + * Electron-phonon and phonon-electron scattering rates by Wannier interpolation * Electron-phonon scattering within the electron-phonon averaged (EPA) approximation - * Polar correction and boundary scattering contributions to transport * Electronic transport coefficients (mobility, conductivity, thermal conductivity, and Seebeck coefficient) #### Phonon Transport - * 3-phonon scattering from thirdOrder.py/ShengBTE or Phono3py force constants - * Boundary and isotope scattering contributions to transport - * Phonon (lattice) thermal conductivity - * Lattice thermal conductivity calculations including ph-el lifetimes + * Phonon (lattice) thermal conductivity, including: + * 3-phonon scattering from thirdOrder.py/ShengBTE or Phono3py force constants + * Boundary, isotope, and phonon-electron scattering contributions + * Lattice thermal conductivity calculations including both ph-ph and ph-el scattering #### And more... - * BTE solutions by RTA, iterative, variational, and relaxons solvers + * BTE solutions by RTA, iterative, variational, and relaxon solvers * Calculation of electron and phonon linewidths or relaxation times on a path - * Wigner transport equation correction for electrons and phonons (Zener tunneling contribution to electron transport) + * Wigner transport equation correction for electrons and phonons * Hydrodynamic transport properties (viscosity) for electrons and phonons diff --git a/doc/sphinx/source/inputFiles.rst b/doc/sphinx/source/inputFiles.rst index 6829d64e..327f7a65 100644 --- a/doc/sphinx/source/inputFiles.rst +++ b/doc/sphinx/source/inputFiles.rst @@ -135,6 +135,8 @@ Phonon BTE Solver * :ref:`dimensionality` +* :ref:`thickness` + * :ref:`constantRelaxationTime` * :ref:`withIsotopeScattering` @@ -222,6 +224,8 @@ Electron BTE Solver * :ref:`dimensionality` +* :ref:`thickness` + * :ref:`constantRelaxationTime` * :ref:`convergenceThresholdBTE` @@ -295,6 +299,10 @@ EPA Transport * :ref:`epaEnergyRange` +* :ref:`dimensionality` + +* :ref:`thickness` + * :ref:`kMesh` * :ref:`temperatures` @@ -482,7 +490,7 @@ Electron Lifetimes on a Path ----------------------------------- -Phonon Dos +Phonon DoS ---------- **Functionality:** Compute the phonon density of states. @@ -1114,7 +1122,7 @@ convergenceThresholdBTE dimensionality ^^^^^^^^^^^^^^ -* **Description:** Input the dimensionality of the material. As a result, transport coefficients tensors will be of size (dim x dim), and units will be suitably scaled for the desired dimensionality. +* **Description:** Input the dimensionality of the material. As a result, transport coefficient tensors will be of size (dim x dim), and units will be suitably scaled for the desired dimensionality. For 2D materials, note that Phoebe assumes the material is oriented in x-y, with height in z. * **Format:** *int* @@ -1123,6 +1131,20 @@ dimensionality * **Default:** `3` +.. _thickness: + +thickness +^^^^^^^^^^^^^^ + +* **Description:** Input the thickness of a 2D material. This is only used if :ref:`dimensionality` = 2, and results in a scaling factor of (cell height/thickness) applied to the transport results. + +* **Format:** *double* + +* **Required:** no + +* **Default:** `1.` + + .. _constantRelaxationTime: constantRelaxationTime diff --git a/doc/sphinx/source/introduction.rst b/doc/sphinx/source/introduction.rst index e9fe5615..afa171b7 100644 --- a/doc/sphinx/source/introduction.rst +++ b/doc/sphinx/source/introduction.rst @@ -3,7 +3,7 @@ Introduction Phoebe is a software for the ab-initio prediction of transport properties, such as electron-phonon limited conductivity and phonon thermal conductivity. -Currently, we support integration with the ab-initio software suite Quantum ESPRESSO for electron-phonon properties and anharmonic force constants through ShengBTE. Phoebe also accepts anharmonic force constants from Phono3py, which interfaces with most mainstream DFT packages. Support for more coming soon. +Currently, we support integration with the ab-initio software suite Quantum ESPRESSO for electron-phonon properties and anharmonic force constants through ShengBTE and Phono3py. Additionally, Phoebe is written in C++ and designed for use on modern HPC infrastructures through hybrid MPI/OpenMP parallelization, distributed memory computation via ScaLAPACK, and support for GPU acceleration using Kokkos. @@ -20,17 +20,19 @@ J. Phys. Mater. 5 035003. (2022). * Electron-phonon limited electronic transport coefficients, including the electrical conductivity, Seebeck coefficient, and electronic thermal conductivity -* Lattice thermal conductivity, including effects from 3-phonon scattering, phonon-isotope and phonon-boundary scattering +* Lattice thermal conductivity, including contributions from 3-phonon, phonon-isotope, phonon-electron, and phonon-boundary scattering -* Electron-phonon scattering through Wannier interpolation, or EPA approximation +* BTE solutions by RTA, iterative, variational, and relaxon solvers -* Calculation of electron and phonon lifetimes (linewidths) +* Electron-phonon transport via Wannier interpolation or EPA approximation -* Electron and phonon viscosities +* Wigner transport equation correction for electrons and phonons + +* Calculation of electron and phonon lifetimes/linewidths (including projected onto a band path) -* Plots of electronic band structure and phonon dispersion relations +* Electron and phonon viscosities -* Plots of electron or phonon density of states +* Plots of electron band structures/phonon dispersion, electron and phonon density of states This documentation contains a brief description of the formalism used to compute these quantities. We expect that the user is already familiar with ab-initio codes and with the fundamentals of solid state physics. A good introduction to the topic is the book "Electrons and Phonons: The Theory of Transport Phenomena in Solids" by Ziman. In the theory section, we add references to contemporary literature for the latest method/algorithm developments. diff --git a/doc/sphinx/source/postProcessing.rst b/doc/sphinx/source/postProcessing.rst index 3987a8db..fc520c8f 100644 --- a/doc/sphinx/source/postProcessing.rst +++ b/doc/sphinx/source/postProcessing.rst @@ -85,6 +85,15 @@ Plots lifetimes vs energy. Example usage:: See note above under ``tau_path.py`` for more info about calcIndex (the 0 at the end). +epa_tau.py +^^^^^^^^^^^^^^^^^^^^^^^^^ + +Plots lifetimes vs energy from an EPA calculation. Example usage:: + + python epa_tau.py epa_relaxation_times.json 0 + +See note above under ``tau_path.py`` for more info about calcIndex (the 0 at the end). + tau_bandResolved.py ^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/sphinx/source/theory/phononBTE.rst b/doc/sphinx/source/theory/phononBTE.rst index 41a37d08..b6cbee68 100644 --- a/doc/sphinx/source/theory/phononBTE.rst +++ b/doc/sphinx/source/theory/phononBTE.rst @@ -437,5 +437,5 @@ Hence, one can work with the same techniques detailed above, provided that we wo * Disadvantage 2: note that the symmetric matrix gains two Cartesian coordinate indices. As a result, in the limiting case of no symmetries in the system (only the identity), the matrix :math:`A^{ij}_{\nu^*\nu'^*}` will still be computed on the same number of wavevectors of :math:`A_{\nu\nu'}`, but occupies 3x3 times more memory without adding any information. Therefore, for low-symmetry systems, consider disabling symmetries. -* Disadvantage 3: The symmetries of the BTE are so far not applicable to the variational and relaxons solveres. This is not so much a problem with implementaiton, but instead is because of a need for a derivation of symmetries for these cases. +* Disadvantage 3: The symmetries of the BTE are so far not applicable to the variational and relaxons solvers. This is not so much a problem with implementation, but instead is because of a need for a derivation of symmetries for these cases. diff --git a/doc/sphinx/source/tutorials/elEpaTransport.rst b/doc/sphinx/source/tutorials/elEpaTransport.rst index 1481a50e..881cb2b6 100644 --- a/doc/sphinx/source/tutorials/elEpaTransport.rst +++ b/doc/sphinx/source/tutorials/elEpaTransport.rst @@ -229,7 +229,7 @@ To do this, let's have a look at the input file ``qeToPhoebeEPA.in``:: electronFourierCutoff = 4. epaMinEnergy = -4. eV epaMaxEnergy = 10. eV - epaNumBins = 10 + epaNumBins = 40 epaSmearingEnergy = 0.05 eV The parameters in this input file are as follows: @@ -278,6 +278,7 @@ Below is an example input file for computing electronic transport properties:: appName = "transportEpa" + useSymmetries = true electronH0Name = "qe-elph/out/silicon.xml", epaFileName = "qe-elph/silicon.phoebe.epa.dat" @@ -293,6 +294,8 @@ The parameters used here are: * :ref:`appName` = `"transportEPA"`: selects the app for computing electronic transport properties with EPA. +* :ref:`useSymmetries`: whether or not to use symmetries of the crystal structure in the calculation. + * :ref:`electronH0Name`: points to the Quantum-ESPRESSO ``*.xml`` file created by ``pw.x``, which contains the electronic single-particle energies. * :ref:`epaFileName`: is the path to the file created at the previous step with ``elPhQeToPhoebe``. @@ -340,16 +343,14 @@ You can learn more about how to post-process these files at :ref:`postprocessing **Files which are always output for this calculation:** -* ``electron_bands.json``: contains the electron band energies used in the calculation. - -* ``electron_dos.json``: contains the electron density of states used in the calculation. +* ``epa_dos.json``: contains the electron density of states used in the calculation. * ``epa_onsager_coefficients.json``: contains the electronic transport coefficients from EPA. * ``epa_relaxation_times.json``: contains the EPA relaxation times at each energy bin value. -To understand how to parse these files in more detail, take a look at the scripts described by the :ref:`postprocessing` page. In particular, if you want to plot lifetimes vs. energy, look at ``tau.py``. If you want to plot the transport coefficients vs. doping or temperature, check out ``transport_coefficients.py``. +To understand how to parse these files in more detail, take a look at the scripts described by the :ref:`postprocessing` page. In particular, if you want to plot lifetimes vs. energy, look at ``epa_tau.py``. If you want to plot the transport coefficients vs. doping or temperature, check out ``transport_coefficients.py``. .. note:: It's a good idea to also use bands.py to plot the band structure. Fourier interpolation of the band structure can be a source of error -- you may need to go to higher k-point meshes to get a reasonable interpolation of the band structure. Check the bands at this stage to make sure they are similar to the true DFT bands -- otherwise, your results could be problematic! diff --git a/example/Silicon-el/reference/variational_onsager_coefficients.json b/example/Silicon-el/reference/variational_onsager_coefficients.json index b5f66c4c..f39df470 100644 --- a/example/Silicon-el/reference/variational_onsager_coefficients.json +++ b/example/Silicon-el/reference/variational_onsager_coefficients.json @@ -10,19 +10,19 @@ "electricalConductivity": [ [ [ - 1067476.083423818, - -289451.1726509079, - 6988.620289292474 + 1067524.755964012, + -289797.97319140926, + 4576.259106139662 ], [ - -288402.6295156577, - 1153574.4448552558, - -339252.03872597613 + -290534.25409045833, + 1153700.6483462334, + -337794.79923207546 ], [ - 6298.688866869923, - -338748.84561951895, - 999624.5190918705 + 2423.3550264169166, + -337602.28847748507, + 999823.4331092425 ] ] ], @@ -30,19 +30,19 @@ "electronicThermalConductivity": [ [ [ - 1.6784141552068075, - -0.17095390966293136, - -0.1925551757224847 + 1.6759588003842663, + -0.1695184449742141, + -0.19316637320993976 ], [ - -0.16887322607943037, - 0.967469430606753, - 0.45201855869035085 + -0.17307030791637715, + 0.9671708308136298, + 0.4527993530482056 ], [ - -0.191703686958177, - 0.4539177794991126, - 1.2702570365076697 + -0.1895368529229464, + 0.4451924672444988, + 1.265051833070335 ] ] ], @@ -50,19 +50,19 @@ "mobility": [ [ [ - 66.6266227276008, - -18.066122864709598, - 0.43619540955680475 + 66.62966062889252, + -18.087768453901255, + 0.28562765358693365 ], [ - -18.00067794376874, - 72.00046026235661, - -21.1744487251345 + -18.133723497245306, + 72.00833726542098, + -21.083494981541044 ], [ - 0.3931332732677861, - -21.14304187885133, - 62.3916607941002 + 0.15125393775778945, + -21.071479404221538, + 62.404076031683935 ] ] ], @@ -71,19 +71,19 @@ "seebeckCoefficient": [ [ [ - -3.834412383472706, - -24.797489282120214, - 14.089965245893998 + -3.7355359725395547, + -24.81297218056885, + 14.075023704624254 ], [ - -14.256947447646306, - -19.694668984349708, - -14.015560582347534 + -14.206137023556979, + -19.708154573007306, + -14.119749174882436 ], [ - 14.752677089503676, - -8.432328471018659, - -59.49649274673085 + 15.116531799699265, + -8.807694880433802, + -59.56690767248755 ] ] ], diff --git a/example/Silicon-epa/reference/epa_dos.json b/example/Silicon-epa/reference/epa_dos.json new file mode 100644 index 00000000..459e5351 --- /dev/null +++ b/example/Silicon-epa/reference/epa_dos.json @@ -0,0 +1,1209 @@ +{ + "dos": [ + 0.0022852187985158193, + 0.0022944645662989526, + 0.002304431107938969, + 0.0023151184234358616, + 0.0023265265127896304, + 0.002338655376000273, + 0.002350297941762473, + 0.0023617650334609203, + 0.002374703390341741, + 0.0023881444809992773, + 0.002403651274076728, + 0.002422372820346082, + 0.0024443091198073402, + 0.0024694601724605, + 0.0024947661527920206, + 0.002521196853850212, + 0.002549475618916622, + 0.002579602447991256, + 0.0026115773410741107, + 0.002645399834237183, + 0.002677115389047235, + 0.0027029375796312504, + 0.0027228664059892062, + 0.002711358051840728, + 0.0026788973839612946, + 0.0026510037238428137, + 0.002627677071485284, + 0.0026089174268886985, + 0.002594724790053067, + 0.0025850991609783773, + 0.0025800405396646403, + 0.0025795489261118488, + 0.002585008644529757, + 0.0026048755283338957, + 0.0026176737765032606, + 0.00261982570815665, + 0.0026129616536627095, + 0.0026150349055375326, + 0.0026186902021535056, + 0.0026239275435106135, + 0.0026307469296088628, + 0.002639148360448266, + 0.002649131836028812, + 0.0026606973563504987, + 0.0026738449214133244, + 0.0026885745312173054, + 0.0027048861857624238, + 0.0027227798850486835, + 0.002742174997677715, + 0.0027623102940358643, + 0.0027831904342717284, + 0.0028413186147160747, + 0.002951240011384291, + 0.003112954624276375, + 0.0032717160265478555, + 0.003097855300903517, + 0.003036655624739618, + 0.0029938308050765387, + 0.0029605286356714266, + 0.002929748937714925, + 0.002901388253901695, + 0.002871200739671667, + 0.002843084611291179, + 0.0028170398687602386, + 0.0027930665120788354, + 0.002771164541246969, + 0.002751333956264656, + 0.0027335747571318766, + 0.0027178869438486417, + 0.0027042705164149485, + 0.002693549056582468, + 0.0026849573939040205, + 0.002676316102038496, + 0.002667625180985913, + 0.0026591778648282828, + 0.0026506214648855733, + 0.0026418750124110064, + 0.0026329385074045646, + 0.0026238119498662644, + 0.0026124686986629462, + 0.0025997190282619127, + 0.002586801625480707, + 0.002573716490319352, + 0.002560486471783936, + 0.002547535186710039, + 0.0025350220166190265, + 0.0025229479128716266, + 0.00251150009746337, + 0.002500818662560505, + 0.002490903608163023, + 0.0024840470176722786, + 0.0024782713489199634, + 0.002473345595494054, + 0.002469269757394554, + 0.0024660438346214603, + 0.0024636678271747763, + 0.0024631079393587453, + 0.002463369167600098, + 0.002464122697498321, + 0.0024653727176681636, + 0.0024671274429224122, + 0.002469387089820521, + 0.0024721516583624905, + 0.0026968927727768323, + 0.0024418739392613395, + 0.0025151312329526445, + 0.003034429045683461, + 0.0032152025872888537, + 0.0025647501389955986, + 0.0025597010821964694, + 0.0025546393601950947, + 0.0025495649729914737, + 0.002544477920585593, + 0.00253937820297747, + 0.002534265820167096, + 0.0025291407721544697, + 0.002524003058939591, + 0.0025188526805224646, + 0.002513689636903095, + 0.0025085139280814676, + 0.0025033255540575926, + 0.002498124514831472, + 0.0024929108104030998, + 0.00248827729261504, + 0.002483973919932755, + 0.0024799760339384177, + 0.002476283634632031, + 0.002472896722013609, + 0.0027608966879738316, + 0.002694660708038876, + 0.0026415409321882224, + 0.0026015373604218917, + 0.0025746499927398665, + 0.0025606610402572174, + 0.0025504135497229656, + 0.0025393292474187633, + 0.002527408310551504, + 0.002514650739121187, + 0.0025010565331278155, + 0.002486625692571384, + 0.002471358217451899, + 0.0024552541077693542, + 0.0024383133635237546, + 0.0024205621093508507, + 0.0024022955218468635, + 0.0023893583533669822, + 0.0023792574810939287, + 0.002368849804355709, + 0.0023581353231523316, + 0.00234711403748379, + 0.002335785947350087, + 0.00232415105275122, + 0.0023122093536871906, + 0.0022999608501579963, + 0.002287405542163645, + 0.002307010672329682, + 0.0024680400459258403, + 0.002500029396669648, + 0.0025131373888541805, + 0.0025098600209145262, + 0.0024904240212312236, + 0.0024548293898042782, + 0.00240307612663368, + 0.002335164231719437, + 0.0022510937050615387, + 0.00218261086680854, + 0.0021283579006754778, + 0.002074762597290411, + 0.0020218249566533505, + 0.0019695449787642877, + 0.0019179226636232257, + 0.0018669580112301634, + 0.0018166510215850999, + 0.0017670016946880387, + 0.0017177891528156704, + 0.0016693751141056315, + 0.0016221007281618264, + 0.0015759659949842516, + 0.0015399820333053965, + 0.0015142610374232585, + 0.0014893332044295274, + 0.0014651985343242066, + 0.0014418570271072953, + 0.001419308682778794, + 0.0013975535013387018, + 0.0013765914827870188, + 0.001356361096353644, + 0.0013369926902499236, + 0.0013185304622149662, + 0.0013009744122487723, + 0.0012843245403513447, + 0.001268580846522679, + 0.0012975021705813616, + 0.0013509784289295664, + 0.0013631099187373682, + 0.0013343742539813751, + 0.001263583837837386, + 0.0011800614832204995, + 0.0011690332832369969, + 0.0011582871026884739, + 0.0011477375117759923, + 0.0011372011939329652, + 0.0011266719931334263, + 0.001116149909377374, + 0.0011056349426648074, + 0.001095127092995729, + 0.0010846263603701382, + 0.0010741327447880353, + 0.0010636462462494175, + 0.0010531668647542874, + 0.001042694600302644, + 0.001032229452894489, + 0.0010217714225298187, + 0.0010113205092086384, + 0.0010008767129309436, + 0.000990450494397441, + 0.0009807305753297448, + 0.0009708845184686043, + 0.00095987594012854, + 0.0009485610855238768, + 0.0009369399546546153, + 0.0009250125475207562, + 0.000912778864122297, + 0.0009002389044592414, + 0.0008873926685315876, + 0.0008742401563393335, + 0.0008607813678824835, + 0.0008469931062870625, + 0.0008327042659508696, + 0.0008181886870318184, + 0.0008043288738450051, + 0.0007904832863147082, + 0.0007766519244409344, + 0.0007628346228451985, + 0.0007488105781775275, + 0.0007343820280079945, + 0.0007195489723366022, + 0.0007043114111633469, + 0.0006886693444882326, + 0.0006726227723112559, + 0.0006563915166275816, + 0.0006407875516207435, + 0.0006239712521894989, + 0.0006059426183338428, + 0.0005867016500537801, + 0.0005662483473493053, + 0.0005445827102204248, + 0.0005217047386671314, + 0.0004976144326894322, + 0.0004732688769938498, + 0.00045841621615233557, + 0.000445335492798058, + 0.00043079695848002815, + 0.00041833207799081904, + 0.0004056683878720973, + 0.00039280588812385986, + 0.0003797445787461102, + 0.0003664844597388442, + 0.0003530255311020664, + 0.00033936779283577244, + 0.00032551124493996616, + 0.00031145588741464365, + 0.0002972017202598093, + 0.0002827487434754586, + 0.0002680969570615959, + 0.00025647754111265654, + 0.00024462596224315185, + 0.0002321864596774886, + 0.00021915903341567037, + 0.0002055436834576936, + 0.00019134040980356206, + 0.00017654921245327152, + 0.00016603849952398183, + 0.0001606144916083271, + 0.00015496959531709534, + 0.0001491038106502848, + 0.00014301713760789728, + 0.0001367095761899309, + 0.00013018112639638754, + 0.00012343178822726535, + 0.0001164615616825663, + 0.00010927044676228819, + 0.00010185844346643331, + 9.422555179499941e-05, + 8.637177174798877e-05, + 7.8297103325399e-05, + 7.000154652723258e-05, + 6.148510135348696e-05, + 5.2214244691523044e-05, + 4.3528321188963454e-05, + 3.5632027095727484e-05, + 2.8525362411815102e-05, + 2.2208327137226303e-05, + 1.6680921271961116e-05, + 1.1943144816019518e-05, + 7.994997769401513e-06, + 4.836480132107106e-06, + 2.467591904136296e-06, + 8.883330854890805e-07, + 9.87036761654613e-08, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 1.4927441146062155e-09, + 1.475404204176711e-08, + 4.173701438103581e-08, + 8.244166113241235e-08, + 1.368679822958967e-07, + 2.0501597787148886e-07, + 2.8688564785918886e-07, + 3.8247699225899664e-07, + 4.917900110709124e-07, + 6.148247042949358e-07, + 7.515810719310673e-07, + 9.020591139793059e-07, + 1.066258830439653e-06, + 1.2441802213121075e-06, + 1.4358232865966702e-06, + 1.6411880262933407e-06, + 1.8602744404021188e-06, + 2.0930825289230047e-06, + 2.3396122918559986e-06, + 2.5998637292011008e-06, + 2.87383684095831e-06, + 3.1615316271276256e-06, + 3.462948087709053e-06, + 3.7780862227025847e-06, + 4.106946032108224e-06, + 4.449527515925976e-06, + 4.80583067415583e-06, + 5.175855506797796e-06, + 6.898678723273927e-06, + 1.2702624589416238e-05, + 2.266769772081461e-05, + 3.679389811747027e-05, + 5.508122577937959e-05, + 7.752968070654431e-05, + 0.00010413926289896878, + 0.00013490997235665004, + 0.0001698418090795823, + 0.0002053924016572728, + 0.00023421948507472385, + 0.0002587814575673882, + 0.0002815028119545206, + 0.0003023835482361255, + 0.00032142366641220536, + 0.00033862316648275724, + 0.0003539820484477787, + 0.0003675003123072724, + 0.0003791779580612401, + 0.00038901498570967976, + 0.0003970113952525897, + 0.00040316718668997236, + 0.0004074823600218278, + 0.0004100398210686938, + 0.00041222601017040176, + 0.0004145349471984687, + 0.0004169666321528957, + 0.00041952106503368185, + 0.00042219824584082706, + 0.00042499817457433096, + 0.0004279208512341943, + 0.00043096627582041747, + 0.00043413444833299955, + 0.00043742536877194063, + 0.00044083903713724103, + 0.00044437545342890114, + 0.0004486459920390718, + 0.00045582014188948704, + 0.00046315817253302533, + 0.0004706600839696883, + 0.00047832587619947396, + 0.0004861555492223821, + 0.000494149103038412, + 0.0005023065376475665, + 0.0005106278530498443, + 0.0005191130492452436, + 0.0005277621262337658, + 0.0005365750840154124, + 0.000545551922590182, + 0.0005546926419580736, + 0.0005639972421190875, + 0.0005734657230732263, + 0.0005830980848204876, + 0.0005907361341545066, + 0.00059562331484914, + 0.0006006312996751095, + 0.000605760088632414, + 0.0006110096817210531, + 0.0006163800789410274, + 0.0006218712802923382, + 0.0006274832857749843, + 0.0006332160953889643, + 0.0006390697091342791, + 0.0006450441270109305, + 0.0006511393490189182, + 0.0006573553751582388, + 0.0006636922054288957, + 0.000670149839830888, + 0.0006767282783642157, + 0.0006834275210288776, + 0.0006902475678248755, + 0.000697188418752209, + 0.0007042500738108786, + 0.0007114325330008811, + 0.0007187357963222196, + 0.0007261658318627886, + 0.0007337356900074617, + 0.0007414458301096203, + 0.000749613234551422, + 0.0007594979730423018, + 0.0007712303124898097, + 0.0007848102528939412, + 0.0008002377942547005, + 0.0008172186431395001, + 0.0008352555280006684, + 0.0008561618093655726, + 0.0008713193940308998, + 0.0008798071225115211, + 0.0008884294585822039, + 0.0008971864022429486, + 0.0009059750046282275, + 0.0009141256712800104, + 0.0009224215634930734, + 0.0009308626812674196, + 0.0009394490246030453, + 0.000948180593499954, + 0.0009572385785214624, + 0.0009680251401091478, + 0.0009808357415289248, + 0.0009956703827808029, + 0.0010125290638647792, + 0.001031411784780851, + 0.0010522612882900638, + 0.0010741513362849862, + 0.0010967578377279806, + 0.0011200807926190473, + 0.0011441202009581878, + 0.0011688760627454047, + 0.0011943483779806952, + 0.001220537146664055, + 0.0012474423687954899, + 0.0012750640443749998, + 0.0013034021734025867, + 0.0013324567558782416, + 0.0013622277918019718, + 0.0013927152811737792, + 0.0014239192239936577, + 0.0014558396202616099, + 0.001488476469977632, + 0.001521829773141735, + 0.0015566180322114903, + 0.0015925385913177815, + 0.001629147703421363, + 0.001666445368522241, + 0.0017044315866204072, + 0.001743106357715859, + 0.0017824696818085986, + 0.0020108736775489924, + 0.002330491036578267, + 0.0024049559819080287, + 0.0022342685135383444, + 0.002178388157834294, + 0.0022234810563476373, + 0.0022569660973659933, + 0.002278299227638483, + 0.0023114066629331183, + 0.0023557593102666884, + 0.002398049456966815, + 0.00243827051436438, + 0.002476051802318442, + 0.0025116909578129713, + 0.002545187980847964, + 0.002576542871423429, + 0.002605755629539366, + 0.0026328262551957747, + 0.002658219469521461, + 0.0026825910690414265, + 0.0027059457267023983, + 0.0027282834425043684, + 0.0027498654008054025, + 0.002772625944449351, + 0.002798703533485995, + 0.0028258090583025103, + 0.002853942518898897, + 0.002883103915275148, + 0.0029132932474312764, + 0.0029445105153672675, + 0.0029767557190831296, + 0.0030101401829124978, + 0.003045130887745353, + 0.003158827665968762, + 0.0032721546179216945, + 0.003377773845697189, + 0.0034760013660938703, + 0.0035673180519141946, + 0.003651728953865285, + 0.0037292340719471398, + 0.0037867853466327438, + 0.0038392104793577445, + 0.0038872928190616753, + 0.0039216788620399365, + 0.0039319934340305415, + 0.003936845768003889, + 0.003938717313452741, + 0.003936912753410502, + 0.0039305382904905, + 0.0039190259513776915, + 0.0039024344884113255, + 0.00388229445054001, + 0.003859019227195484, + 0.003832608818377754, + 0.0038030632240868122, + 0.0037517204484383603, + 0.0036886762033664735, + 0.0037816586806912124, + 0.0036529839393109305, + 0.003450769730622205, + 0.0034025947424494085, + 0.003355268861909193, + 0.0033084292223448117, + 0.0032620758237562604, + 0.0032162086661435218, + 0.003170827749506602, + 0.0031259330738455035, + 0.0030815246391602334, + 0.003037602445450787, + 0.0029941664927171543, + 0.0029512167809593462, + 0.002908753310177369, + 0.0028667760803712064, + 0.0028243942103777105, + 0.002781461171685018, + 0.002738013065663962, + 0.002693855365000556, + 0.0026486165471348777, + 0.0026025503179402548, + 0.0025556774042024373, + 0.0025114158716398782, + 0.0024731740302192397, + 0.002435035405649141, + 0.0023969999979295865, + 0.002359067807060578, + 0.0023212388330421034, + 0.002283513075874172, + 0.0022458905355567937, + 0.002208371212089952, + 0.00217095510547365, + 0.0021336422157078845, + 0.0021009259368146415, + 0.0020737474417929116, + 0.002046060233892902, + 0.002017864313114623, + 0.0019874493301375857, + 0.0019564296735234492, + 0.0019206007749233043, + 0.0018851640764936588, + 0.0018498227159103229, + 0.0018135686124131447, + 0.0017859958467512595, + 0.0017590029351182424, + 0.0017325522982485722, + 0.0017066439361422432, + 0.0016812778487992539, + 0.0016564540362196086, + 0.0016321724984033055, + 0.0016084332353503471, + 0.0015852362470607261, + 0.0015625815335344483, + 0.0015404690947715147, + 0.0015188989307719227, + 0.001497871041535669, + 0.0014788774838943013, + 0.0014616935686200933, + 0.0014449998747255016, + 0.0014287964022105257, + 0.0014130831510751698 + ], + "dosUnit": "1/eV", + "energies": [ + 3.2426301970281006, + 3.2526301970281004, + 3.2626301970281006, + 3.2726301970281004, + 3.2826301970281007, + 3.2926301970281004, + 3.3026301970281007, + 3.3126301970281005, + 3.3226301970281007, + 3.3326301970281005, + 3.3426301970281007, + 3.3526301970281005, + 3.3626301970281003, + 3.3726301970281005, + 3.3826301970281003, + 3.3926301970281005, + 3.4026301970281003, + 3.4126301970281, + 3.4226301970281003, + 3.4326301970281, + 3.4426301970281004, + 3.4526301970281, + 3.4626301970281004, + 3.4726301970281, + 3.4826301970281004, + 3.4926301970281, + 3.5026301970281004, + 3.5126301970281, + 3.5226301970281004, + 3.5326301970281, + 3.5426301970281004, + 3.5526301970281002, + 3.5626301970281, + 3.5726301970281003, + 3.5826301970281, + 3.5926301970281003, + 3.6026301970281, + 3.6126301970281003, + 3.6226301970281, + 3.6326301970281003, + 3.6426301970281, + 3.6526301970281003, + 3.6626301970281, + 3.6726301970281003, + 3.6826301970281, + 3.6926301970281004, + 3.7026301970281, + 3.7126301970281004, + 3.7226301970281, + 3.7326301970281004, + 3.7426301970281, + 3.7526301970281004, + 3.7626301970281, + 3.7726301970281004, + 3.7826301970281, + 3.7926301970281004, + 3.8026301970281002, + 3.8126301970281005, + 3.8226301970281003, + 3.8326301970281005, + 3.8426301970281003, + 3.8526301970281005, + 3.8626301970281003, + 3.8726301970281005, + 3.8826301970281003, + 3.8926301970281005, + 3.9026301970281003, + 3.9126301970281006, + 3.9226301970281003, + 3.9326301970281006, + 3.9426301970281004, + 3.9526301970281006, + 3.9626301970281004, + 3.9726301970281006, + 3.9826301970281004, + 3.9926301970281006, + 4.0026301970281, + 4.0126301970281, + 4.022630197028101, + 4.032630197028101, + 4.0426301970281004, + 4.0526301970281, + 4.062630197028101, + 4.072630197028101, + 4.0826301970281005, + 4.0926301970281, + 4.102630197028101, + 4.112630197028101, + 4.1226301970281005, + 4.1326301970281, + 4.142630197028101, + 4.152630197028101, + 4.1626301970281006, + 4.1726301970281, + 4.182630197028101, + 4.192630197028101, + 4.202630197028101, + 4.2126301970281, + 4.222630197028101, + 4.232630197028101, + 4.242630197028101, + 4.2526301970281, + 4.2626301970281, + 4.272630197028101, + 4.282630197028101, + 4.2926301970281004, + 4.3026301970281, + 4.312630197028101, + 4.322630197028101, + 4.3326301970281005, + 4.3426301970281, + 4.352630197028101, + 4.362630197028101, + 4.3726301970281, + 4.3826301970281, + 4.3926301970281, + 4.402630197028101, + 4.4126301970281, + 4.4226301970281, + 4.4326301970281, + 4.442630197028101, + 4.4526301970281, + 4.4626301970281, + 4.4726301970281, + 4.482630197028101, + 4.4926301970281, + 4.5026301970281, + 4.5126301970281, + 4.522630197028101, + 4.5326301970281, + 4.5426301970281004, + 4.5526301970281, + 4.562630197028101, + 4.5726301970281, + 4.5826301970281005, + 4.5926301970281, + 4.602630197028101, + 4.6126301970281, + 4.6226301970281005, + 4.6326301970281, + 4.642630197028101, + 4.6526301970281, + 4.6626301970281006, + 4.6726301970281, + 4.6826301970281, + 4.6926301970281, + 4.7026301970281, + 4.7126301970281, + 4.7226301970281, + 4.7326301970281, + 4.7426301970281, + 4.7526301970281, + 4.7626301970281, + 4.7726301970281, + 4.7826301970281, + 4.7926301970281004, + 4.8026301970281, + 4.8126301970281, + 4.8226301970281, + 4.8326301970281005, + 4.8426301970281, + 4.8526301970281, + 4.8626301970281, + 4.8726301970281005, + 4.8826301970281, + 4.8926301970281, + 4.9026301970281, + 4.9126301970281006, + 4.9226301970281, + 4.9326301970281, + 4.9426301970281, + 4.952630197028101, + 4.9626301970281, + 4.9726301970281, + 4.9826301970281, + 4.992630197028101, + 5.0026301970281, + 5.0126301970281, + 5.0226301970281, + 5.032630197028101, + 5.0426301970281004, + 5.0526301970281, + 5.0626301970281, + 5.072630197028101, + 5.0826301970281005, + 5.0926301970281, + 5.1026301970281, + 5.112630197028101, + 5.1226301970281005, + 5.1326301970281, + 5.1426301970281, + 5.152630197028101, + 5.1626301970281006, + 5.1726301970281, + 5.1826301970281, + 5.192630197028101, + 5.202630197028101, + 5.2126301970281, + 5.2226301970281, + 5.232630197028101, + 5.242630197028101, + 5.2526301970281, + 5.2626301970281, + 5.272630197028101, + 5.282630197028101, + 5.2926301970281004, + 5.3026301970281, + 5.312630197028101, + 5.322630197028101, + 5.3326301970281005, + 5.3426301970281, + 5.3526301970281, + 5.362630197028101, + 5.3726301970281005, + 5.3826301970281, + 5.3926301970281, + 5.402630197028101, + 5.4126301970281006, + 5.4226301970281, + 5.4326301970281, + 5.442630197028101, + 5.452630197028101, + 5.4626301970281, + 5.4726301970281, + 5.4826301970281, + 5.492630197028101, + 5.5026301970280995, + 5.5126301970281, + 5.5226301970281, + 5.532630197028101, + 5.5426301970281, + 5.5526301970281, + 5.5626301970281, + 5.572630197028101, + 5.5826301970281, + 5.5926301970281, + 5.6026301970281, + 5.612630197028101, + 5.6226301970281, + 5.6326301970281, + 5.6426301970281, + 5.652630197028101, + 5.6626301970281, + 5.6726301970281, + 5.6826301970281, + 5.692630197028101, + 5.7026301970281, + 5.7126301970281, + 5.7226301970281, + 5.732630197028101, + 5.7426301970281, + 5.7526301970281, + 5.7626301970281, + 5.772630197028101, + 5.7826301970281, + 5.7926301970281004, + 5.8026301970281, + 5.812630197028101, + 5.8226301970281, + 5.8326301970281005, + 5.8426301970281, + 5.852630197028101, + 5.8626301970281, + 5.8726301970281005, + 5.8826301970281, + 5.892630197028101, + 5.9026301970281, + 5.9126301970281006, + 5.9226301970281, + 5.932630197028101, + 5.9426301970281, + 5.952630197028101, + 5.9626301970281, + 5.972630197028101, + 5.9826301970281, + 5.992630197028101, + 6.0026301970281, + 6.012630197028101, + 6.0226301970281, + 6.032630197028101, + 6.0426301970281004, + 6.052630197028101, + 6.0626301970281, + 6.072630197028101, + 6.0826301970281005, + 6.092630197028101, + 6.1026301970281, + 6.112630197028101, + 6.1226301970281005, + 6.1326301970281, + 6.1426301970281, + 6.1526301970281, + 6.1626301970281006, + 6.1726301970281, + 6.1826301970281, + 6.1926301970281, + 6.202630197028101, + 6.2126301970281, + 6.2226301970281, + 6.2326301970281, + 6.242630197028101, + 6.2526301970281, + 6.2626301970281, + 6.2726301970281, + 6.282630197028101, + 6.2926301970281004, + 6.3026301970281, + 6.3126301970281, + 6.322630197028101, + 6.3326301970281005, + 6.3426301970281, + 6.3526301970281, + 6.362630197028101, + 6.3726301970281005, + 6.3826301970281, + 6.3926301970281, + 6.402630197028101, + 6.4126301970281006, + 6.4226301970281, + 6.4326301970281, + 6.4426301970281, + 6.452630197028101, + 6.4626301970281, + 6.4726301970281, + 6.4826301970281, + 6.492630197028101, + 6.5026301970281, + 6.5126301970281, + 6.5226301970281, + 6.532630197028101, + 6.5426301970281004, + 6.5526301970281, + 6.5626301970281, + 6.572630197028101, + 6.5826301970281005, + 6.5926301970281, + 6.6026301970281, + 6.612630197028101, + 6.6226301970281005, + 6.6326301970281, + 6.6426301970281, + 6.652630197028101, + 6.6626301970281006, + 6.6726301970281, + 6.6826301970281, + 6.692630197028101, + 6.702630197028101, + 6.7126301970281, + 6.7226301970281, + 6.732630197028101, + 6.742630197028101, + 6.7526301970281, + 6.7626301970281, + 6.772630197028101, + 6.782630197028101, + 6.7926301970281004, + 6.8026301970281, + 6.812630197028101, + 6.822630197028101, + 6.8326301970281, + 6.8426301970281, + 6.852630197028101, + 6.862630197028101, + 6.8726301970281, + 6.8826301970281, + 6.892630197028101, + 6.902630197028101, + 6.9126301970281, + 6.9226301970281, + 6.932630197028101, + 6.942630197028101, + 6.9526301970281, + 6.9626301970281, + 6.972630197028101, + 6.982630197028101, + 6.9926301970281, + 7.0026301970281, + 7.012630197028101, + 7.022630197028101, + 7.0326301970281, + 7.0426301970281004, + 7.052630197028101, + 7.062630197028101, + 7.0726301970281, + 7.0826301970281, + 7.0926301970281, + 7.102630197028101, + 7.1126301970281, + 7.1226301970281, + 7.1326301970281, + 7.142630197028101, + 7.1526301970281, + 7.1626301970281, + 7.1726301970281, + 7.182630197028101, + 7.1926301970281, + 7.2026301970281, + 7.2126301970281, + 7.222630197028101, + 7.2326301970281, + 7.242630197028099, + 7.2526301970281, + 7.262630197028101, + 7.2726301970281, + 7.282630197028099, + 7.2926301970281004, + 7.302630197028101, + 7.3126301970281, + 7.322630197028099, + 7.3326301970281005, + 7.342630197028101, + 7.3526301970281, + 7.362630197028099, + 7.3726301970281005, + 7.382630197028101, + 7.3926301970281, + 7.402630197028099, + 7.4126301970281006, + 7.422630197028101, + 7.4326301970281, + 7.442630197028099, + 7.452630197028101, + 7.462630197028101, + 7.4726301970281, + 7.482630197028099, + 7.492630197028101, + 7.502630197028101, + 7.5126301970281, + 7.522630197028099, + 7.5326301970281, + 7.542630197028101, + 7.5526301970281, + 7.562630197028099, + 7.5726301970281, + 7.582630197028101, + 7.5926301970281, + 7.602630197028099, + 7.6126301970281, + 7.622630197028101, + 7.6326301970281, + 7.642630197028099, + 7.6526301970281, + 7.662630197028101, + 7.6726301970281, + 7.682630197028099, + 7.6926301970281, + 7.7026301970281015, + 7.7126301970281, + 7.722630197028099, + 7.7326301970281, + 7.7426301970281015, + 7.7526301970281, + 7.762630197028099, + 7.7726301970281, + 7.782630197028101, + 7.7926301970281004, + 7.802630197028099, + 7.8126301970281, + 7.822630197028101, + 7.8326301970281005, + 7.842630197028099, + 7.8526301970281, + 7.862630197028101, + 7.8726301970281005, + 7.882630197028099, + 7.8926301970281, + 7.902630197028101, + 7.9126301970281006, + 7.9226301970280995, + 7.9326301970281, + 7.942630197028101, + 7.952630197028101, + 7.9626301970280995, + 7.9726301970281, + 7.982630197028101, + 7.992630197028101, + 8.0026301970281, + 8.0126301970281, + 8.022630197028102, + 8.0326301970281, + 8.0426301970281, + 8.052630197028101, + 8.062630197028101, + 8.0726301970281, + 8.0826301970281, + 8.0926301970281, + 8.102630197028102, + 8.1126301970281, + 8.1226301970281, + 8.132630197028101, + 8.142630197028101, + 8.1526301970281, + 8.1626301970281, + 8.1726301970281, + 8.182630197028102, + 8.1926301970281, + 8.2026301970281, + 8.212630197028101, + 8.222630197028101, + 8.2326301970281, + 8.2426301970281, + 8.2526301970281, + 8.262630197028102, + 8.2726301970281, + 8.2826301970281, + 8.292630197028101, + 8.302630197028101, + 8.312630197028101, + 8.3226301970281, + 8.3326301970281, + 8.342630197028102, + 8.3526301970281, + 8.3626301970281, + 8.3726301970281, + 8.3826301970281, + 8.392630197028101, + 8.4026301970281, + 8.412630197028099, + 8.4226301970281, + 8.4326301970281, + 8.4426301970281, + 8.4526301970281, + 8.4626301970281, + 8.472630197028101, + 8.4826301970281, + 8.492630197028099, + 8.5026301970281, + 8.5126301970281, + 8.5226301970281, + 8.5326301970281, + 8.5426301970281, + 8.552630197028101, + 8.562630197028101, + 8.572630197028099, + 8.5826301970281, + 8.5926301970281, + 8.6026301970281, + 8.6126301970281, + 8.6226301970281, + 8.632630197028101, + 8.6426301970281, + 8.652630197028099, + 8.6626301970281, + 8.6726301970281, + 8.6826301970281, + 8.6926301970281, + 8.7026301970281, + 8.712630197028101, + 8.7226301970281, + 8.732630197028099, + 8.7426301970281, + 8.7526301970281, + 8.7626301970281, + 8.7726301970281, + 8.7826301970281, + 8.792630197028101, + 8.8026301970281, + 8.8126301970281, + 8.8226301970281, + 8.8326301970281, + 8.8426301970281, + 8.8526301970281, + 8.8626301970281, + 8.872630197028101, + 8.8826301970281, + 8.8926301970281, + 8.9026301970281, + 8.9126301970281, + 8.9226301970281, + 8.9326301970281, + 8.9426301970281, + 8.952630197028101, + 8.9626301970281, + 8.9726301970281, + 8.9826301970281, + 8.9926301970281, + 9.0026301970281, + 9.0126301970281, + 9.0226301970281, + 9.032630197028102, + 9.0426301970281, + 9.0526301970281, + 9.062630197028101, + 9.0726301970281, + 9.0826301970281, + 9.0926301970281, + 9.1026301970281, + 9.112630197028102, + 9.1226301970281, + 9.1326301970281, + 9.142630197028101, + 9.1526301970281, + 9.1626301970281, + 9.1726301970281, + 9.1826301970281, + 9.192630197028102, + 9.2026301970281, + 9.2126301970281, + 9.222630197028101, + 9.2326301970281 + ], + "energyUnit": "eV", + "particleType": "electron" +} diff --git a/example/Silicon-ph/reference/variational_phonon_thermal_cond.json b/example/Silicon-ph/reference/variational_phonon_thermal_cond.json index b8f4c5f3..a70bc0d2 100644 --- a/example/Silicon-ph/reference/variational_phonon_thermal_cond.json +++ b/example/Silicon-ph/reference/variational_phonon_thermal_cond.json @@ -8,18 +8,18 @@ [ [ 60.97660348006533, - 0.05939673195170657, - 0.05510061800860491 + 0.05939673195170365, + 0.055100618008605204 ], [ - 0.030536633906653364, + 0.030536633906649856, 60.698641062099384, - -0.002692541767271135 + -0.0026925417672721587 ], [ - 0.043433199458606485, - -0.024502698690385013, - 60.73017036224241 + 0.04343319945860528, + -0.024502698690383698, + 60.730170362242426 ] ] ], diff --git a/scripts/plotScripts/epa_tau.py b/scripts/plotScripts/epa_tau.py new file mode 100644 index 00000000..4b270983 --- /dev/null +++ b/scripts/plotScripts/epa_tau.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +import json +import matplotlib.pyplot as plt +import numpy as np +import argparse +import itertools +import os +import sys + +# script to plot the EPA relaxation times found in epa_relaxation_times.json + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="Plot relaxation times that " + "have been generated by Phoebe using the EPA method.") + parser.add_argument("INPUT", + help="Name of the EPA JSON file with relaxation times") + parser.add_argument("calcIndex", + help="Index of temperature/chemical potential to plot", + default=0) + args = parser.parse_args() + + # we select one calculation and check that the json file is correct + try: + calcIndex = int(args.calcIndex) + except ValueError: + raise ValueError("calcIndex should be an integer") + + # load in the json output + jfileName = args.INPUT + with open(jfileName) as jfile: + data = json.load(jfile) + + try: + data['relaxationTimes'] + except KeyError: + raise KeyError("relaxation times not found." + "Are you using the correct input json file?") + + particleType = data['particleType'] + + # unpack the json file + tau = np.array(data['relaxationTimes']) # dimensions (iCalc, ik, ib) + tau[np.where(tau==None)] = 0 # remove None values (from gamma pt acoustic ph) + energies = np.array(data['energies']) # dimensions (iCalc, ik, ib) + linewidths = np.array(data['linewidths']) # dimensions (iCalc, ik, ib) + mu = np.array(data['chemicalPotentials']) + T = np.array(data['temperatures']) + + linewidths = linewidths[calcIndex].flatten() + energies = energies[calcIndex].flatten() + tau = tau[calcIndex].flatten() + mu = mu[calcIndex] + energies = energies - mu + + print("Calculation Temperature: ", T[calcIndex]) + + for y, name in [[tau,'tau'], [linewidths,'Gamma']]: + + # plot the lifetimes + plt.figure(figsize=(5,5)) + + plt.scatter(energies, y, marker='o', s=18, color='royalblue') + + # plot aesthetics + plt.yscale('log') + plt.xlabel(r'Energy [' + data['energyUnit'] +']',fontsize=12) + if name == 'tau': + units = ' [' + data['relaxationTimeUnit'] + ']' + else: + units = ' [' + data['linewidthsUnit'] + ']' + + plt.ylabel(r'$\{}_{{'.format(name) + data['particleType'] + '}$' + + units, fontsize=12) + + if (particleType=="phonon"): + plt.xlim(0,None) + + # Find limits of the y axis + zeroIndex = np.argwhere(y<=0.) + y = np.delete(y, zeroIndex) + ymin = 10**np.floor(np.log10(np.min(y))) + ymax = 10**np.ceil(np.log10(np.max(y))) + plt.ylim(ymin, ymax) + + plt.tight_layout() + plt.ylim(0.1, 1000) + plotFileName = os.path.splitext(jfileName)[0] + ".{}.png".format(name.lower()) + plt.savefig(plotFileName,dpi=150) diff --git a/scripts/plotScripts/readCouplingElPh.py b/scripts/plotScripts/readCouplingElPh.py new file mode 100644 index 00000000..1040a132 --- /dev/null +++ b/scripts/plotScripts/readCouplingElPh.py @@ -0,0 +1,27 @@ +import numpy as np +import h5py + +# read in the file +f = h5py.File('coupling.elph.phoebe.hdf5', 'r') + +# k and q are reconstructed from this below +points = np.array(f["pointsPairsCrystal"]) +points = np.reshape(points,(6,-1)).T +numPointsPairs = points.shape[0] + +# unfold the points, which are listed in pairs as kx ky kz qx qy qz +ks = points[:,0:3] +qs = points[:,3:6] + +# set up band info +bandRange1 = np.array(f["elBandRange1"]) +bandRange2 = np.array(f["elBandRange2"]) +bandRangePh = np.array(f["phModeRange"]) +numBands1 = bandRange1[1] - bandRange1[0] + 1 +numBands2 = bandRange2[1] - bandRange2[0] + 1 +numBandsPh = bandRangePh[1] - bandRangePh[0] + 1 + +couplingMat = np.array(f["elphCouplingMat"]) +couplingMat= np.reshape(couplingMat, (numPointsPairs, numBands1, numBands2, numBandsPh)) + +print(couplingMat) diff --git a/scripts/plotScripts/tau.py b/scripts/plotScripts/tau.py index 43901ddd..478efbf4 100755 --- a/scripts/plotScripts/tau.py +++ b/scripts/plotScripts/tau.py @@ -7,40 +7,24 @@ import os import sys -def noWindow(data,jfileName,calcIndex): - # unpack the json file - tau = np.array(data['relaxationTimes']) # dimensions (iCalc, ik, ib) - tau[np.where(tau==None)] = 0 # remove None values (from gamma pt acoustic ph) - energies = np.array(data['energies']) # dimensions (iCalc, ik, ib) - linewidths = np.array(data['linewidths']) # dimensions (iCalc, ik, ib) - nbands = energies.shape[2] - mu = np.array(data['chemicalPotentials']) - T = np.array(data['temperatures']) - - linewidths = linewidths[calcIndex].flatten() - energies = energies[calcIndex].flatten() - tau = tau[calcIndex].flatten() - mu = mu[calcIndex] +def parseData(data,jfileName,calcIndex): - print("Calculation Temperature: ", T[calcIndex]) - - return energies-mu, tau, linewidths - -def withWindow(data,jfileName,calcIndex): # unpack the json file - srcTau = np.array(data['relaxationTimes']) # dimensions (iCalc, ik, ib) - srcTau[np.where(srcTau==None)] = 0 # remove None values (from gamma pt acoustic ph) - srcEnergies = np.array(data['energies']) # dimensions (iCalc, ik, ib) - srcLinewidths = np.array(data['linewidths']) # dimensions (iCalc, ik, ib) tau = [] energies = [] linewidths = [] - for x in srcTau: + + for x in data['relaxationTimes']: tau.append(np.array(list(itertools.chain(*x)))) - for x in srcEnergies: + tau = np.array(tau) + tau[np.where(tau==None)] = 0 # remove None values due to acoustic modes + + for x in data['energies']: energies.append(np.array(list(itertools.chain(*x)))) - for x in srcLinewidths: + + for x in data['linewidths']: linewidths.append(np.array(list(itertools.chain(*x)))) + mu = np.array(data['chemicalPotentials']) T = np.array(data['temperatures']) @@ -52,11 +36,8 @@ def withWindow(data,jfileName,calcIndex): parser = argparse.ArgumentParser(description="Plot relaxation times that " "have been generated by Phoebe") - parser.add_argument("INPUT", - help="Name of the JSON file with relaxation times") - parser.add_argument("calcIndex", - help="Index of temperature/chemical potential to plot", - default=0) + parser.add_argument("INPUT", help="Name of the JSON file with relaxation times") + parser.add_argument("calcIndex", help="Index of temperature/chemical potential to plot", default=0) args = parser.parse_args() # we select one calculation and check that the json file is correct @@ -78,20 +59,14 @@ def withWindow(data,jfileName,calcIndex): particleType = data['particleType'] - # if a window was used to filter points in the calculation, - # the format of the output can be slightly different - if len(np.array(data['relaxationTimes']).shape) == 3: - energies, tau, linewidths = noWindow(data,jfileName,calcIndex) - else: - energies, tau, linewidths = withWindow(data,jfileName,calcIndex) - + # parse data from json files + energies, tau, linewidths = parseData(data, jfileName, calcIndex) for y, name in [[tau,'tau'], [linewidths,'Gamma']]: - # plot the lifetimes plt.figure(figsize=(5,5)) - plt.scatter(energies, y, marker='o', s=18, color='royalblue') + plt.scatter(energies, y, marker='o', s=18, color='yellowgreen') # plot aesthetics plt.yscale('log') @@ -101,8 +76,7 @@ def withWindow(data,jfileName,calcIndex): else: units = ' [' + data['linewidthsUnit'] + ']' - plt.ylabel(r'$\{}_{{'.format(name) + data['particleType'] + '}$' + - units, fontsize=12) + plt.ylabel(r'$\{}_{{'.format(name) + data['particleType'] + '}$' + units, fontsize=12) if (particleType=="phonon"): plt.xlim(0,None) @@ -117,3 +91,5 @@ def withWindow(data,jfileName,calcIndex): plt.tight_layout() plotFileName = os.path.splitext(jfileName)[0] + ".{}.pdf".format(name.lower()) plt.savefig(plotFileName) + + diff --git a/src/apps/bands_app.cpp b/src/apps/bands_app.cpp index e736c801..9afca800 100644 --- a/src/apps/bands_app.cpp +++ b/src/apps/bands_app.cpp @@ -41,6 +41,7 @@ void PhononBandsApp::run(Context &context) { if (mpi->mpiHead()) { std::cout << "Finishing phonon bands calculation" << std::endl; } + } /* --------------------- ElectronWannierBandsApp --------------------- */ @@ -154,9 +155,9 @@ void outputBandsToJSON(FullBandStructure &fullBandStructure, Context &context, auto pathLabels = context.getPathLabels(); // check if this point is one of the high sym points, // and if it is, save the index - if (coord[0] == extremaCoordinates[extremaCount][0] && - coord[1] == extremaCoordinates[extremaCount][1] && - coord[2] == extremaCoordinates[extremaCount][2]) { + if ( abs(coord[0] - extremaCoordinates[extremaCount][0]) < 1e-8 && + abs(coord[1] - extremaCoordinates[extremaCount][1]) < 1e-8 && + abs(coord[2] - extremaCoordinates[extremaCount][2]) < 1e-8 ) { pathLabelIndices.push_back(ik); // if the next extrema is the same as this one, add the index diff --git a/src/apps/dos_app.h b/src/apps/dos_app.h index 35238b3b..b6ee63e8 100644 --- a/src/apps/dos_app.h +++ b/src/apps/dos_app.h @@ -22,4 +22,14 @@ class ElectronFourierDosApp : public App { void checkRequirements(Context &context) override; }; +/** + * Helper function to output dos to a json file + * @param energies: list of energies to write to file + * @param dos: list of dos values to write to file + * @param particle: electron or phonon particle object + * @param outFileName: name of the dos .json file + */ +void outputDOSToJSON(std::vector energies, std::vector dos, + Particle& particle, const std::string &outFileName); + #endif diff --git a/src/apps/electron_wannier_transport_app.cpp b/src/apps/electron_wannier_transport_app.cpp index 3501f29d..1c78109a 100644 --- a/src/apps/electron_wannier_transport_app.cpp +++ b/src/apps/electron_wannier_transport_app.cpp @@ -431,8 +431,7 @@ void ElectronWannierTransportApp::runVariationalMethod( auto outF = scatteringMatrix.dot(inF); VectorBTE az2E = outF[0]; VectorBTE az2T = outF[1]; - transportCoefficients.calcVariational(az2E, az2T, zNewE, zNewT, bE, bT, - preconditioning); + transportCoefficients.calcVariational(az2E, az2T, zNewE, zNewT, bE, bT, preconditioning); transportCoefficients.print(iter); elCond = transportCoefficients.getElectricalConductivity(); thCond = transportCoefficients.getThermalConductivity(); @@ -440,7 +439,7 @@ void ElectronWannierTransportApp::runVariationalMethod( // decide whether to exit or run the next iteration double deltaE = findMaxRelativeDifference(elCond, elCondOld); double deltaT = findMaxRelativeDifference(thCond, thCondOld); - if ((deltaE < threshold) && (deltaT < threshold)) { + if ((deltaE < threshold) && (deltaT < threshold) && iter > 2) { // recompute our final guess for transport coefficients transportCoefficients.calcFromSymmetricPopulation(zNewE, zNewT); break; diff --git a/src/apps/phonon_transport_app.cpp b/src/apps/phonon_transport_app.cpp index b7a0b5b5..89cbb823 100644 --- a/src/apps/phonon_transport_app.cpp +++ b/src/apps/phonon_transport_app.cpp @@ -251,6 +251,9 @@ void PhononTransportApp::run(Context &context) { std::cout << "Starting variational BTE solver\n" << std::endl; } + // NOTE this is currently un-preconditioned! We should update the algorithm as here + // See numerical recipies: the minimal residual algorithm, 2.7 Sparse Linear Systems, pg 89 + // note: each iteration should take approximately twice as long as // the iterative method above (in the way it's written here. diff --git a/src/apps/transport_epa_app.cpp b/src/apps/transport_epa_app.cpp index 05692f46..11b358e3 100644 --- a/src/apps/transport_epa_app.cpp +++ b/src/apps/transport_epa_app.cpp @@ -14,6 +14,7 @@ #include "utilities.h" #include #include +#include "dos_app.h" void TransportEpaApp::run(Context &context) { @@ -30,10 +31,10 @@ void TransportEpaApp::run(Context &context) { // in principle, we should add 1 to account for ends of energy interval // I will not do that, because will work with the centers of energy steps - int numEnergies = int((maxEnergy - minEnergy) / energyStep); + size_t numEnergies = size_t((maxEnergy - minEnergy) / energyStep); // energies at the centers of energy steps Eigen::VectorXd energies(numEnergies); - for (int i = 0; i < numEnergies; ++i) { + for (size_t i = 0; i < numEnergies; ++i) { // add 0.5 to be in the middle of the energy step energies(i) = (double(i) + 0.5) * energyStep + minEnergy; } @@ -42,9 +43,7 @@ void TransportEpaApp::run(Context &context) { Points fullPoints(crystal, context.getKMesh()); fullPoints.setIrreduciblePoints(); - if (mpi->mpiHead()) { - std::cout << "\nBuilding electronic band structure" << std::endl; - } + // Building electronic band structure ----------------------- // filter to only the bands relevant to transport electronH0.trimBands(context, minEnergy, maxEnergy); @@ -60,7 +59,10 @@ void TransportEpaApp::run(Context &context) { if (mpi->mpiHead()) { std::cout << "Starting EPA with " << numEnergies << " energies and " - << bandStructure.getNumStates() << " states" << std::endl; + << bandStructure.getNumStates() << " states." << std::endl; + std::cout << std::setprecision(4) << "Energy grid ranges from " << minEnergy*energyRyToEv << " to " + << maxEnergy*energyRyToEv + << " eV with " << numEnergies << " bins." << std::endl; } //-------------------------------- @@ -68,11 +70,19 @@ void TransportEpaApp::run(Context &context) { auto t2 = calcEnergyProjVelocity(context, bandStructure, energies); Eigen::Tensor energyProjVelocity = std::get<0>(t2); Eigen::VectorXd dos = std::get<1>(t2); - + { + // output to file (convert to std::vector for the sake of reusing dos output function from dos_app) + // braces to make this go out of scope after writing + std::vector dosVec(dos.data(), dos.data() + dos.rows() * dos.cols()); + std::vector energyVec(energies.data(), energies.data() + energies.rows() * energies.cols()); + Particle particle = bandStructure.getParticle(); + outputDOSToJSON(energyVec, dosVec, particle, "epa_dos.json"); + } //-------------------------------- // Calculate EPA scattering rates VectorEPA scatteringRates = getScatteringRates(context, statisticsSweep, energies, crystal, dos); + outputToJSON("epa_relaxation_times.json", scatteringRates, statisticsSweep, energies, context); @@ -183,7 +193,7 @@ VectorEPA TransportEpaApp::getScatteringRates( return crtRate; } - int spinFactor = 2; + double spinFactor = 2.; if (context.getHasSpinOrbit()) { spinFactor = 1; } @@ -212,7 +222,7 @@ VectorEPA TransportEpaApp::getScatteringRates( } } - LoopPrint loopPrint("calculation of EPA scattering rates", "phonon modes", + LoopPrint loopPrint("calculation of EPA scattering rates", "energy bins per MPI process", int(mpi->divideWorkIter(numEnergies).size())); VectorEPA epaRate(statisticsSweep, numEnergies, 1); diff --git a/src/bands/bandstructure.cpp b/src/bands/bandstructure.cpp index 570a07a9..01f65de3 100644 --- a/src/bands/bandstructure.cpp +++ b/src/bands/bandstructure.cpp @@ -325,7 +325,7 @@ void FullBandStructure::setEnergies(Point &point, Eigen::VectorXd &energies_) { int ik = point.getIndex(); if (!energies.indicesAreLocal(0,ik)) { // col distributed, only need to check ik - Error("Cannot access a non-local energy"); + Error("DeveloperError: Cannot access a non-local energy in setEnergies."); } for (int ib = 0; ib < energies.localRows(); ib++) { energies(ib, ik) = energies_(ib); @@ -334,9 +334,11 @@ void FullBandStructure::setEnergies(Point &point, Eigen::VectorXd &energies_) { void FullBandStructure::setVelocities( Point &point, Eigen::Tensor, 3> &velocities_) { + if (!hasVelocities) { - Error("FullBandStructure was initialized without velocities"); + Error("FullBandStructure was initialized without velocities, cannot set velocities."); } + // we convert from a tensor to a vector (how it's stored in memory) Eigen::VectorXcd tmpVelocities_(numBands * numBands * 3); for (int i = 0; i < numBands; i++) { @@ -349,10 +351,12 @@ void FullBandStructure::setVelocities( } } int ik = point.getIndex(); + if(ik >= numPoints || ik < 0) std::cout << "found a bad ik value " << ik << std::endl; if (!velocities.indicesAreLocal(0,ik)) { // col distributed, only need to check ik - Error("Cannot access a non-local velocity"); + Error("DeveloperError: Cannot set a non-local velocity in distributed velocity vector."); } + // here this isn't a band index, it's actually an index over all compressed band indices for (int ib = 0; ib < velocities.localRows(); ib++) { velocities(ib, ik) = tmpVelocities_(ib); } diff --git a/src/bte/el_scattering.cpp b/src/bte/el_scattering.cpp index bfb8f890..1f5b4687 100644 --- a/src/bte/el_scattering.cpp +++ b/src/bte/el_scattering.cpp @@ -109,7 +109,8 @@ void ElScatteringMatrix::builder(VectorBTE *linewidth, mpi->allReduceSum(&innerFermi); if (smearing->getType() == DeltaFunction::tetrahedron) { - Error("Tetrahedron method not supported by electron scattering"); + Error("Tetrahedron method not supported by electron scattering. It will not work properly with a\n" + "filtering window, and we will almost always need to use this for electrons."); // that's because it doesn't work with the window the way it's implemented, // and we will almost always have a window for electrons } diff --git a/src/bte/scattering.cpp b/src/bte/scattering.cpp index 8dfe2076..9429dc59 100644 --- a/src/bte/scattering.cpp +++ b/src/bte/scattering.cpp @@ -38,10 +38,13 @@ ScatteringMatrix::ScatteringMatrix(Context &context_, excludeIndices.push_back(iBte); } - Eigen::Vector3d k = outerBandStructure.getWavevector(isIdx); - if (k.squaredNorm() > 1e-8 && en < 0.) { - Warning("Found a phonon mode q!=0 with negative energies. " - "Consider improving the quality of your DFT phonon calculation.\n"); + Eigen::Vector3d q = outerBandStructure.getWavevector(isIdx); + if (q.squaredNorm() > 1e-8 && en < 0.) { + q = outerBandStructure.getPoints().cartesianToCrystal(q); + Warning("Found a phonon mode q!=0 with negative energy.\n" + "Consider improving the quality of your DFT phonon calculation, or applying a different sum rule.\n" + "Energy: " + std::to_string(en * energyRyToEv*100.) + " meV, q (in crystal) = " + + std::to_string(q(0)) + ", " + std::to_string(q(1)) + ", " + std::to_string(q(2))); } } } @@ -69,10 +72,12 @@ ScatteringMatrix::ScatteringMatrix(Context &context_, } smearing = DeltaFunction::smearingFactory(context, innerBandStructure); - if ( // innerBandStructure != outerBandStructure && - smearing->getType() == DeltaFunction::tetrahedron) { - Error("Tetrahedron smearing for transport untested and thus blocked"); - // not for linewidths. Although this should be double-checked + // block for electron scattering matrix + if ( (innerBandStructure.getParticle().isElectron() + || outerBandStructure.getParticle().isElectron() ) + && smearing->getType() == DeltaFunction::tetrahedron) { + Error("Tetrahedron smearing can cause issues for smearing with a state filtering window," + "\nwhich one essentially will always want to use for electrons."); } } diff --git a/src/common_kokkos.cpp b/src/common_kokkos.cpp index 989f99aa..1a9715e3 100644 --- a/src/common_kokkos.cpp +++ b/src/common_kokkos.cpp @@ -102,23 +102,13 @@ void kokkosZHEEV(StridedComplexView3D &A, DoubleView2D &W) { DeviceManager *kokkosDeviceMemory = nullptr; DeviceManager::DeviceManager() { + memoryUsed = 0.; char *memStr = std::getenv("MAXMEM"); if (memStr != nullptr) { memoryTotal = std::atof(memStr) * 1.0e9; } else { memoryTotal = 16.0e9; // 16 Gb is our educated guess for available memory - -// size_t freeMemory, totalMemory; -//#if defined(KOKKOS_ENABLE_SERIAL) || defined(KOKKOS_ENABLE_OPENMP) -// -//#elif defined(KOKKOS_ENABLE_CUDA) -// cudaMemGetInfo (&freeMemory, &totalMemory ); -// memoryUsed += totalMemory - freeMemory; -// memoryTotal = double(totalMemory); -//#else -// Error("Implement DevicecManager for this backend"); -//#endif } } diff --git a/src/context.cpp b/src/context.cpp index df3aae3d..db05c9b5 100644 --- a/src/context.cpp +++ b/src/context.cpp @@ -83,6 +83,12 @@ double parseDoubleWithUnits(std::string &line) { if (patternInString(line, "mum")) { x /= distanceBohrToMum; } + if (patternInString(line, "ang")) { + x /= distanceBohrToAng; + } + if (patternInString(line, "Bohr") || patternInString(line, "bohr")) { + x /= 1.; // it's already in atomic units + } return x; } @@ -133,7 +139,6 @@ std::tuple parseDoubleVectorComponent(std::string line) { return std::make_tuple(idx,val); } - /** Parse a string of format "key = value units" to return an integer value. */ int parseInt(std::string &line) { @@ -548,6 +553,9 @@ void Context::setupFromInput(const std::string &fileName) { if (parameterName == "dimensionality") { dimensionality = parseInt(val); } + if (parameterName == "thickness") { + thickness = parseDoubleWithUnits(val); + } if (parameterName == "dosMinEnergy") { dosMinEnergy = parseDoubleWithUnits(val); } @@ -760,6 +768,40 @@ void Context::setupFromInput(const std::string &fileName) { lineCounter += 1; } } + +void Context::inputSanityCheck() { + + // disallow symmetries when relaxons are used + for (const std::string &s : solverBTE) { + if (s.compare("variational") == 0 || s.compare("relaxons") == 0) { + if(useSymmetries) { + Error("Variational and relaxons solvers cannot be used with symmetries!"); + } + } + + // Warn users about relaxons with even meshes + if (s.compare("relaxons") == 0) { + if((qMesh.prod() % 2 == 0 && qMesh.prod() != 0) || (kMesh.prod() % 2 == 0 && kMesh.prod() != 0)) { + Warning("Relaxons solver should be run with an odd points mesh for reasons of eigenvector parity."); + } + if(smearingMethod == 1) { // adaptive smearing + Warning("Adaptive smearing can cause problems with the quality of the scattering matrix, and may cause negative eigenvalues to appear.\n" + "If this occurs, you should switch to Gaussian."); + } + } + } + + // warn the user if thickness is not set but 2d is + if(dimensionality == 2) { + if(thickness == 1.) { + Warning("You have set dimensionality = 2 but not the thickness input parameter.\n" + "This means your final result will have (height of cell / 1.) applied, when it should be\n" + "height of cell / thickness!"); + } + } + +} + // helper functions for printInputSummary template void printVector(const std::string& varName, std::vector vec) { @@ -790,6 +832,7 @@ void Context::printInputSummary(const std::string &fileName) { // crystal structure parameters ------------------- std::cout << "useSymmetries = " << useSymmetries << std::endl; std::cout << "dimensionality = " << dimensionality << std::endl; + if(dimensionality != 3) std::cout << "thickness = " << thickness * distanceBohrToAng << " ang" << std::endl; std::cout << std::endl; // phonon parameters ------------------------------- @@ -974,8 +1017,7 @@ void Context::printInputSummary(const std::string &fileName) { if (dopings.size() != 0) printVectorXd("dopings", dopings, "cm^-3"); if (chemicalPotentials.size() != 0) - printVectorXd("chemicalPotentials", chemicalPotentials * energyRyToEv, - "eV"); + printVectorXd("chemicalPotentials", chemicalPotentials * energyRyToEv, "eV"); if (!std::isnan(minChemicalPotential)) std::cout << "minChemicalPotential = " << minChemicalPotential * energyRyToEv << " eV" << std::endl; @@ -989,12 +1031,13 @@ void Context::printInputSummary(const std::string &fileName) { std::cout << "eFermiRange = " << eFermiRange << " eV" << std::endl; if (!std::isnan(fermiLevel)) std::cout << "fermiLevel = " << fermiLevel * energyRyToEv << std::endl; - if (!std::isnan(numOccupiedStates)) + if (!std::isnan(numOccupiedStates)) { if(!hasSpinOrbit) { // need to account for spin factor std::cout << "numOccupiedStates = " << numOccupiedStates/2.0 << std::endl; } else { std::cout << "numOccupiedStates = " << numOccupiedStates << std::endl; } + } } // should not be printed when phellifetimes app is run if (appName.find("honon") != std::string::npos && appName.find("lectron") == std::string::npos) { @@ -1112,6 +1155,9 @@ void Context::printInputSummary(const std::string &fileName) { std::cout << "numOccupiedStates = " << numOccupiedStates << std::endl; std::cout << "---------------------------------------------\n" << std::endl; } + + inputSanityCheck(); // call the input sanity check to make sure these make sense + } std::string Context::getPhFC2FileName() { return phFC2FileName; } @@ -1202,6 +1248,8 @@ int Context::getMaxIterationsBTE() const { return maxIterationsBTE; } int Context::getDimensionality() const { return dimensionality; } +double Context::getThickness() const { return thickness; } + double Context::getDosMinEnergy() const { return dosMinEnergy; } double Context::getDosMaxEnergy() const { return dosMaxEnergy; } diff --git a/src/context.h b/src/context.h index a8eb95d8..45000e14 100644 --- a/src/context.h +++ b/src/context.h @@ -52,6 +52,7 @@ class Context { bool hasSpinOrbit = false; int dimensionality = 3; + double thickness = 1.; // material thickness or cross area for lower dimensions double dosMinEnergy = std::numeric_limits::quiet_NaN(); double dosMaxEnergy = std::numeric_limits::quiet_NaN(); @@ -75,8 +76,7 @@ class Context { Eigen::VectorXd customIsotopeCouplings; Eigen::VectorXd customMasses; - // add RTA boundary scattering in phonon scattering matrix - // boundary length for isotope scattering + // add RTA boundary scattering in scattering matrix double boundaryLength = std::numeric_limits::quiet_NaN(); std::string elphFileName; @@ -266,6 +266,8 @@ class Context { int getDimensionality() const; + double getThickness() const; + double getDosMinEnergy() const; double getDosMaxEnergy() const; @@ -345,6 +347,11 @@ class Context { */ void printInputSummary(const std::string &fileName); + /** Sanity checks the input variables to make sure they agree. + * TODO: This function should be replaced with a better system in the future. + */ + void inputSanityCheck(); + Eigen::VectorXi getCoreElectrons(); void setCoreElectrons(const Eigen::VectorXi &x); diff --git a/src/crystal.cpp b/src/crystal.cpp index a8e52525..bae4b66b 100644 --- a/src/crystal.cpp +++ b/src/crystal.cpp @@ -637,3 +637,13 @@ Crystal::buildWignerSeitzVectorsWithShift(const Eigen::Vector3i &grid, } return std::make_tuple(bravaisVectors, degeneracies); } + +// change of basis methods +Eigen::Vector3d Crystal::crystalToCartesian(const Eigen::Vector3d &point) { + return getReciprocalUnitCell() * point; +} + +Eigen::Vector3d Crystal::cartesianToCrystal(const Eigen::Vector3d &point) { + Eigen::Vector3d p = getReciprocalUnitCell().inverse() * point; + return p; +} \ No newline at end of file diff --git a/src/crystal.h b/src/crystal.h index fdd84a02..7be67abe 100644 --- a/src/crystal.h +++ b/src/crystal.h @@ -171,6 +171,18 @@ class Crystal { const Eigen::VectorXd &getAtomicIsotopeCouplings(); + /** Change of basis method to convert crystal wavevector to cartesian coordinates. + * @param point: a point in crystal coordinates + * @return the point in cartesian coordiantes + * */ + Eigen::Vector3d crystalToCartesian(const Eigen::Vector3d &point); + + /** Change of basis method to convert crystal wavevector to cartesian coordinates. + * @param point: a point in cartesian coordinates + * @return the point in crystal coordiantes + * */ + Eigen::Vector3d cartesianToCrystal(const Eigen::Vector3d &point); + /** Build the list of Bravais lattice vectors (real space) that live within * the Wigner Seitz zone of a super cell * which is grid(0) x grid(1) x grid(2) bigger than the unitCell. diff --git a/src/harmonic/electron_h0_fourier.cpp b/src/harmonic/electron_h0_fourier.cpp index b1e47456..790dac15 100644 --- a/src/harmonic/electron_h0_fourier.cpp +++ b/src/harmonic/electron_h0_fourier.cpp @@ -47,6 +47,7 @@ ElectronH0Fourier::ElectronH0Fourier(Crystal &crystal_, const Points& coarsePoin double cutoff_) : crystal(crystal_), coarseBandStructure(coarseBandStructure_), coarsePoints(coarsePoints_), particle(Particle::electron) { + numBands = coarseBandStructure.getNumBands(); cutoff = cutoff_; numDataPoints = coarseBandStructure.getNumPoints(); @@ -73,39 +74,6 @@ ElectronH0Fourier::ElectronH0Fourier(Crystal &crystal_, const Points& coarsePoin loopPrint.close(); } -// Copy constructor -ElectronH0Fourier::ElectronH0Fourier(const ElectronH0Fourier &that) - : crystal(that.crystal), coarseBandStructure(that.coarseBandStructure), - coarsePoints(that.coarsePoints), particle(that.particle), - expansionCoefficients(that.expansionCoefficients), - numBands(that.numBands), cutoff(that.cutoff), - numDataPoints(that.numDataPoints), - numPositionVectors(that.numPositionVectors), - minDistance(that.minDistance), - positionDegeneracies(that.positionDegeneracies), - positionVectors(that.positionVectors), refWavevector(that.refWavevector) { -} - -// Copy assignment -ElectronH0Fourier &ElectronH0Fourier::operator=(const ElectronH0Fourier &that) { - if (this != &that) { - crystal = that.crystal; - coarseBandStructure = that.coarseBandStructure; - coarsePoints = that.coarsePoints; - particle = that.particle; - expansionCoefficients = that.expansionCoefficients; - numBands = that.numBands; - cutoff = that.cutoff; - numDataPoints = that.numDataPoints; - numPositionVectors = that.numPositionVectors; - minDistance = that.minDistance; - positionDegeneracies = that.positionDegeneracies; - positionVectors = that.positionVectors; - refWavevector = that.refWavevector; - } - return *this; -} - Particle ElectronH0Fourier::getParticle() { return particle; } int ElectronH0Fourier::getNumBands() { return numBands; } @@ -156,13 +124,21 @@ FullBandStructure ElectronH0Fourier::populate(Points &fullPoints, const bool &withVelocities, const bool &withEigenvectors, const bool isDistributed) { + FullBandStructure fullBandStructure(numBands, particle, withVelocities, withEigenvectors, fullPoints, isDistributed); - LoopPrint loopPrint("populating electronic band structure", "bands", numBands); + LoopPrint loopPrint("populating electronic band structure", "states", + fullBandStructure.getNumStates()); + + // this is mpi parallel in wavevector indices + std::vector kIndices = fullBandStructure.getWavevectorIndices(); + for (size_t i = 0; i(tup); @@ -446,7 +422,7 @@ ElectronH0Fourier::getGroupVelocityFromCoordinates(Eigen::Vector3d &wavevector, } StridedComplexView3D ElectronH0Fourier::kokkosBatchedBuildBlochHamiltonian( - const DoubleView2D &cartesianCoordinates) { + [[maybe_unused]] const DoubleView2D &cartesianCoordinates) { Error("Kokkos not implemented in ElectronH0Fourier"); // dummy return to silence a warning StridedComplexView3D temp; @@ -455,7 +431,7 @@ StridedComplexView3D ElectronH0Fourier::kokkosBatchedBuildBlochHamiltonian( std::tuple ElectronH0Fourier::kokkosBatchedDiagonalizeWithVelocities( - const DoubleView2D &cartesianCoordinates) { + [[maybe_unused]] const DoubleView2D &cartesianCoordinates) { Error("Kokkos not implemented in ElectronH0Fourier"); std::tuple temp; return temp; @@ -463,7 +439,8 @@ ElectronH0Fourier::kokkosBatchedDiagonalizeWithVelocities( std::tuple ElectronH0Fourier::kokkosBatchedDiagonalizeFromCoordinates( - const DoubleView2D &cartesianCoordinates, const bool withMassScaling) { + [[maybe_unused]] const DoubleView2D &cartesianCoordinates, + [[maybe_unused]] const bool withMassScaling) { Error("Kokkos not implemented in ElectronH0Fourier"); std::tuple temp; return temp; diff --git a/src/harmonic/electron_h0_fourier.h b/src/harmonic/electron_h0_fourier.h index a9decc2f..fbf2ae5d 100644 --- a/src/harmonic/electron_h0_fourier.h +++ b/src/harmonic/electron_h0_fourier.h @@ -32,14 +32,6 @@ class ElectronH0Fourier : public HarmonicHamiltonian { ElectronH0Fourier(Crystal &crystal_, const Points& coarsePoints_, const FullBandStructure& coarseBandStructure_, double cutoff); - /** Copy constructor - */ - ElectronH0Fourier(const ElectronH0Fourier &that); - - /** Copy assignment - */ - ElectronH0Fourier &operator=(const ElectronH0Fourier &that); - /** Method to return that the underlying is that of an electronic Fermion. */ Particle getParticle() override; diff --git a/src/harmonic/phonon_h0.cpp b/src/harmonic/phonon_h0.cpp index f9cb43f0..45e50b27 100644 --- a/src/harmonic/phonon_h0.cpp +++ b/src/harmonic/phonon_h0.cpp @@ -10,11 +10,16 @@ #include "points.h" #include "utilities.h" +#ifdef HDF5_AVAIL +#include +#endif + + PhononH0::PhononH0(Crystal &crystal, const Eigen::Matrix3d &dielectricMatrix_, const Eigen::Tensor &bornCharges_, Eigen::Tensor &forceConstants_, const std::string &sumRule) - : particle(Particle::phonon) { + : particle(Particle::phonon), crystal(crystal) { // in this section, we save as class properties a few variables // that are needed for the diagonalization of phonon frequencies @@ -245,7 +250,7 @@ PhononH0::PhononH0(Crystal &crystal, const Eigen::Matrix3d &dielectricMatrix_, // copy constructor PhononH0::PhononH0(const PhononH0 &that) : particle(that.particle), hasDielectric(that.hasDielectric), - numAtoms(that.numAtoms), numBands(that.numBands), + numAtoms(that.numAtoms), numBands(that.numBands), crystal(that.crystal), volumeUnitCell(that.volumeUnitCell), atomicSpecies(that.atomicSpecies), speciesMasses(that.speciesMasses), atomicPositions(that.atomicPositions), dielectricMatrix(that.dielectricMatrix), bornCharges(that.bornCharges), @@ -274,6 +279,7 @@ PhononH0 &PhononH0::operator=(const PhononH0 &that) { hasDielectric = that.hasDielectric; numAtoms = that.numAtoms; numBands = that.numBands; + crystal = that.crystal; volumeUnitCell = that.volumeUnitCell; atomicSpecies = that.atomicSpecies; speciesMasses = that.speciesMasses; @@ -287,7 +293,6 @@ PhononH0 &PhononH0::operator=(const PhononH0 &that) { mat2R = that.mat2R; gVectors = that.gVectors; longRangeCorrection1 = that.longRangeCorrection1; - atomicMasses_d = that.atomicMasses_d; longRangeCorrection1_d = that.longRangeCorrection1_d; gVectors_d = that.gVectors_d; @@ -331,6 +336,7 @@ PhononH0::diagonalize(Point &point) { return std::make_tuple(energies, eigenvectors); } +// TODO why not just make mass scaling = true the default, then eliminate this function? std::tuple PhononH0::diagonalizeFromCoordinates(Eigen::Vector3d &q) { bool withMassScaling = true; @@ -704,7 +710,6 @@ void PhononH0::reorderDynamicalMatrix(const Eigen::Matrix3d& directUnitCell, } // next, we reorder the dynamical matrix along the bravais lattice vectors - bravaisVectors = Eigen::MatrixXd::Zero(3, numBravaisVectors); weights = Eigen::VectorXd::Zero(numBravaisVectors); mat2R.resize(3, 3, numAtoms, numAtoms, numBravaisVectors); @@ -767,6 +772,7 @@ void PhononH0::reorderDynamicalMatrix(const Eigen::Matrix3d& directUnitCell, void PhononH0::shortRangeTerm(Eigen::Tensor, 4> &dyn, const Eigen::VectorXd &q) { + // calculates the dynamical matrix at q from the (short-range part of the) // force constants21, by doing the Fourier transform of the force constants @@ -1078,3 +1084,70 @@ int PhononH0::getIndexEigenvector(const int &iAt, const int &iPol, return compress2Indices(iAt, iPol, nAtoms, 3); } +void PhononH0::printDynToHDF5(Eigen::Vector3d& qCrys) { + + // wavevector enters in crystal but must be used in cartesian + auto qCart = crystal.crystalToCartesian(qCrys); + + // Construct dynmat for this point + Eigen::Tensor, 4> dyn; + dyn.resize(3, 3, numAtoms, numAtoms); + dyn.setZero(); // might be unnecessary + + // first, the short range term, which is just a Fourier transform + shortRangeTerm(dyn, qCart); + + // then the long range term, which uses some convergence + // tricks by X. Gonze et al. + if (hasDielectric) { + addLongRangeTerm(dyn, qCart); + } + + // the contents can be written to file like this + #ifdef HDF5_AVAIL + + if (mpi->mpiHead()) { + + try { + + std::string outFileName = "dynamical_matrix.hdf5"; + + // if the hdf5 file is there already, we want to delete it. Occasionally + // these files seem to get stuck open when a process dies while writing to + // them, (even if a python script dies) and then they can't be overwritten + // properly. + std::remove(&outFileName[0]); + + // open the hdf5 file + HighFive::File file(outFileName, HighFive::File::Overwrite); + + // flatten the tensor in a vector + Eigen::VectorXcd dynFlat = Eigen::Map(dyn.data(), dyn.size()); + + // write dyn to hdf5 + HighFive::DataSet ddyn = file.createDataSet>( + "/dynamicalMatrix", HighFive::DataSpace::From(dynFlat)); + ddyn.write(dynFlat); + + // write the qpoint + HighFive::DataSet dq = file.createDataSet("/crystalWavevector", HighFive::DataSpace::From(qCrys)); + dq.write(qCrys); + + // write the qpoint + Eigen::VectorXi dims(4); + dims(0) = 3; + dims(1) = 3; + dims(2) = numAtoms; + dims(3) = numAtoms; + HighFive::DataSet ddims = file.createDataSet("/dimensions", HighFive::DataSpace::From(dims)); + ddims.write(dims); + + } catch (std::exception &error) { + Error("Issue writing dynamical matrix to hdf5."); + } + } + #else + Error("One needs to build with HDF5 to write the dynamical matrix to file!"); + #endif + +} diff --git a/src/harmonic/phonon_h0.h b/src/harmonic/phonon_h0.h index 8cd8b3e8..1879dd0c 100644 --- a/src/harmonic/phonon_h0.h +++ b/src/harmonic/phonon_h0.h @@ -147,6 +147,13 @@ class PhononH0 : public HarmonicHamiltonian { * call of the kokkosBatched functions. */ int estimateBatchSize(const bool& withVelocity) override; + + /** Helper function to print the dynamical matrix file to HDF5, for developer testing purposes. + * @param qCrys: a 3d eigen vector in crystal coordinates of the phonon wavevector. + * @param points: the points object with with q belongs to. Here, used only to convert q to cartesian internally. + */ + void printDynToHDF5(Eigen::Vector3d& qCrys); + protected: /** Impose the acoustic sum rule on force constants and Born charges * @param sumRule: name of the sum rule to be used @@ -161,6 +168,7 @@ class PhononH0 : public HarmonicHamiltonian { const Eigen::Tensor& forceConstants); Particle particle; + Crystal &crystal; bool hasDielectric = false; int numAtoms; @@ -216,18 +224,23 @@ class PhononH0 : public HarmonicHamiltonian { /** Adds the long range correction to the dynamical matrix due to dipole-ion * interaction. + * @param dyn: the dynamical matrix in the shape 3,3,natoms,natoms + * @param q: a 3d eigen vector with the cartesian coordinates of the phonon wavevector. */ void addLongRangeTerm(Eigen::Tensor, 4> &dyn, const Eigen::VectorXd &q); /** This part computes the slow-range part of the dynamical matrix, which is * the Fourier transform of the force constants. + * @param dyn: the dynamical matrix in the shape 3,3,natoms,natoms + * @param q: a 3d eigen vector with the cartesian coordinates of the phonon wavevector. */ void shortRangeTerm(Eigen::Tensor, 4> &dyn, const Eigen::VectorXd &q); /** dynDiagonalize diagonalizes the dynamical matrix and returns eigenvalues and * eigenvectors. + * @param dyn: the dynamical matrix in the shape 3,3,natoms,natoms */ std::tuple dynDiagonalize( Eigen::Tensor, 4> &dyn); diff --git a/src/observable/phonon_thermal_cond.cpp b/src/observable/phonon_thermal_cond.cpp index 4b0fdf0f..6c852f61 100644 --- a/src/observable/phonon_thermal_cond.cpp +++ b/src/observable/phonon_thermal_cond.cpp @@ -10,13 +10,26 @@ #include PhononThermalConductivity::PhononThermalConductivity( - Context &context_, StatisticsSweep &statisticsSweep_, Crystal &crystal_, - BaseBandStructure &bandStructure_) - : Observable(context_, statisticsSweep_, crystal_), - bandStructure(bandStructure_) { - tensordxd = - Eigen::Tensor(numCalculations, dimensionality, dimensionality); + Context &context_, StatisticsSweep &statisticsSweep_, Crystal &crystal_, BaseBandStructure &bandStructure_) + : Observable(context_, statisticsSweep_, crystal_), bandStructure(bandStructure_) { + + tensordxd = Eigen::Tensor(numCalculations, dimensionality, dimensionality); tensordxd.setZero(); + + // set up units for writing to file + thCondUnits = "W /(m K)"; + if (dimensionality == 3) { + thCondConversion = thConductivityAuToSi; + } else if (dimensionality == 2) { + // multiply by the height of the cell / thickness of the cell to convert 3D -> 2D. + // Because the unit cell volume is already reduced for dimensionality, + // we only need to divide by thickness. + //double height = crystal.getDirectUnitCell()(2,2); + thCondConversion = thConductivityAuToSi * (1. / context.getThickness()); + } else { // dim = 1 + Warning("1D conductivity should be manually adjusted for the cross-section of the material."); + thCondConversion = thConductivityAuToSi; + } } // copy constructor @@ -50,8 +63,7 @@ void PhononThermalConductivity::calcFromCanonicalPopulation(VectorBTE &f) { void PhononThermalConductivity::calcFromPopulation(VectorBTE &n) { - double norm = 1. / context.getQMesh().prod() / - crystal.getVolumeUnitCell(dimensionality); + double norm = 1. / context.getQMesh().prod() / crystal.getVolumeUnitCell(dimensionality); auto excludeIndices = n.excludeIndices; @@ -60,7 +72,7 @@ void PhononThermalConductivity::calcFromPopulation(VectorBTE &n) { std::vector iss = bandStructure.parallelIrrStateIterator(); int niss = iss.size(); - Kokkos::View> tensordxd_k(tensordxd.data(), numCalculations, 3, 3); + Kokkos::View> tensordxd_k(tensordxd.data(), numCalculations, dimensionality, dimensionality); Kokkos::Experimental::ScatterView scatter_tensordxd(tensordxd_k); Kokkos::parallel_for("phonon_thermal_cond", Kokkos::RangePolicy(0, niss), [&] (int iis){ @@ -86,13 +98,13 @@ void PhononThermalConductivity::calcFromPopulation(VectorBTE &n) { for (int iCalc = 0; iCalc < statisticsSweep.getNumCalculations(); iCalc++) { Eigen::Vector3d nRot; - for (int i : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { nRot(i) = n(iCalc, i, iBte); } nRot = rot * nRot; - for (int j : {0, 1, 2}) { - for (int i : {0, 1, 2}) { + for (int j = 0; j < dimensionality; j++) { + for (int i = 0; i < dimensionality; i++) { tensorPrivate(iCalc, i, j) += nRot(i) * vel(j) * en * norm; } } @@ -135,6 +147,7 @@ void PhononThermalConductivity::calcFromPopulation(VectorBTE &n) { void PhononThermalConductivity::calcVariational(VectorBTE &af, VectorBTE &f, VectorBTE &b) { + double norm = 1. / context.getQMesh().prod() / crystal.getVolumeUnitCell(dimensionality); auto excludeIndices = f.excludeIndices; @@ -149,8 +162,9 @@ void PhononThermalConductivity::calcVariational(VectorBTE &af, VectorBTE &f, std::vector iss = bandStructure.parallelIrrStateIterator(); int niss = iss.size(); - Kokkos::View> y1_k(y1.data(), numCalculations, 3, 3), - y2_k(y2.data(), numCalculations, 3, 3); + // TODO should these be dimensionality? + Kokkos::View> y1_k(y1.data(), numCalculations, dimensionality, dimensionality), + y2_k(y2.data(), numCalculations, dimensionality, dimensionality); Kokkos::Experimental::ScatterView scatter_y1(y1_k), scatter_y2(y2_k); Kokkos::parallel_for("electron_viscosity", Kokkos::RangePolicy(0, niss), [&] (int iis){ auto x1 = scatter_y1.access(); @@ -158,8 +172,7 @@ void PhononThermalConductivity::calcVariational(VectorBTE &af, VectorBTE &f, int is = iss[iis]; // skip the acoustic phonons - if (std::find(excludeIndices.begin(), excludeIndices.end(), is) != - excludeIndices.end()) { + if (std::find(excludeIndices.begin(), excludeIndices.end(), is) != excludeIndices.end()) { return; } @@ -177,7 +190,7 @@ void PhononThermalConductivity::calcVariational(VectorBTE &af, VectorBTE &f, double norm2 = norm * temp * temp; Eigen::Vector3d fRot, afRot, bRot; - for (int i : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { fRot(i) = f(iCalc, i, isBte); afRot(i) = af(iCalc, i, isBte); bRot(i) = b(iCalc, i, isBte); @@ -186,8 +199,8 @@ void PhononThermalConductivity::calcVariational(VectorBTE &af, VectorBTE &f, afRot = rot * afRot; bRot = rot * bRot; - for (int i : {0, 1, 2}) { - for (int j : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { + for (int j = 0; j < dimensionality; j++) { x1(iCalc, i, j) += fRot(i) * afRot(j) * norm2; x2(iCalc, i, j) += fRot(i) * bRot(j) * norm2; } @@ -247,7 +260,7 @@ void PhononThermalConductivity::calcFromRelaxons( int iCalc = 0; // relaxons only allows one calc in memory double temp = statisticsSweep.getCalcStatistics(iCalc).temperature; double chemPot = statisticsSweep.getCalcStatistics(iCalc).chemicalPotential; - VectorBTE population(statisticsSweep, bandStructure, 3); + VectorBTE population(statisticsSweep, bandStructure, dimensionality); // if we only calculated some eigenvalues, // we should not include any alpha @@ -323,7 +336,7 @@ void PhononThermalConductivity::calcFromRelaxons( shared(bandStructure, eigenvectors, relaxonsPopulation, population, \ scatteringMatrix, iCalc, localStates, nlocalStates) { - Eigen::MatrixXd popPrivate(3, bandStructure.irrStateIterator().size()); + Eigen::MatrixXd popPrivate(dimensionality, bandStructure.irrStateIterator().size()); popPrivate.setZero(); #pragma omp for nowait @@ -342,8 +355,8 @@ void PhononThermalConductivity::calcFromRelaxons( #pragma omp critical for (int is = 0; is < popPrivate.cols(); is++) { - for (int iDim = 0; iDim < 3; iDim++) { - population(iCalc, iDim, is) += popPrivate(iDim, is); + for (int i = 0; i < dimensionality; i++) { + population(iCalc, i, is) += popPrivate(i, is); } } } @@ -351,12 +364,12 @@ void PhononThermalConductivity::calcFromRelaxons( } else { // case without symmetries ------------------------------------------ - Eigen::MatrixXd relaxonsPopulation(eigenvalues.size(), 3); + Eigen::MatrixXd relaxonsPopulation(eigenvalues.size(), dimensionality); relaxonsPopulation.setZero(); #pragma omp parallel default(none) \ shared(eigenvectors, eigenvalues, temp, chemPot, particle, relaxonsPopulation, localStates, nlocalStates) { - Eigen::MatrixXd relaxonsPopPrivate(eigenvalues.size(), 3); + Eigen::MatrixXd relaxonsPopPrivate(eigenvalues.size(), dimensionality); relaxonsPopPrivate.setZero(); #pragma omp for nowait for(int ilocalState = 0; ilocalState < nlocalStates; ilocalState++){ @@ -369,7 +382,7 @@ void PhononThermalConductivity::calcFromRelaxons( auto vel = bandStructure.getGroupVelocity(isIdx); double term = sqrt(particle.getPopPopPm1(en, temp, chemPot)); double dndt = particle.getDndt(en, temp, chemPot); - for (int i = 0; i < 3; i++) { + for (int i = 0; i < dimensionality; i++) { relaxonsPopPrivate(alpha, i) += dndt / term * vel(i) * eigenvectors(is, alpha) / eigenvalues(alpha); @@ -378,7 +391,7 @@ void PhononThermalConductivity::calcFromRelaxons( } #pragma omp critical for (int alpha = 0; alpha < eigenvalues.size(); alpha++) { - for (int i = 0; i < 3; i++) { + for (int i = 0; i < dimensionality; i++) { relaxonsPopulation(alpha, i) += relaxonsPopPrivate(alpha, i); } } @@ -389,7 +402,7 @@ void PhononThermalConductivity::calcFromRelaxons( #pragma omp parallel default(none) \ shared(eigenvectors, relaxonsPopulation, population, iCalc, localStates, nlocalStates) { - Eigen::MatrixXd popPrivate(3, bandStructure.getNumStates()); + Eigen::MatrixXd popPrivate(dimensionality, bandStructure.getNumStates()); popPrivate.setZero(); #pragma omp for nowait @@ -397,13 +410,13 @@ void PhononThermalConductivity::calcFromRelaxons( std::tuple tup0 = localStates[ilocalState]; int is = std::get<0>(tup0); int alpha = std::get<1>(tup0); - for (int i = 0; i < 3; i++) { + for (int i = 0; i < dimensionality; i++) { popPrivate(i, is) += eigenvectors(is, alpha) * relaxonsPopulation(alpha, i); } } #pragma omp critical for (int is = 0; is < bandStructure.getNumStates(); is++) { - for (int i = 0; i < 3; i++) { + for (int i = 0; i < dimensionality; i++) { population(iCalc, i, is) += popPrivate(i, is); } } @@ -423,8 +436,8 @@ void PhononThermalConductivity::calcFromRelaxons( int iBte = bandStructure.stateToBte(isIdx).get(); if (en > 0.) { double term = particle.getPopPopPm1(en, temp, chemPot); - for (int iDim = 0; iDim < 3; iDim++) { - population(iCalc, iDim, iBte) *= sqrt(term); + for (int i = 0; i < dimensionality; i++) { + population(iCalc, i, iBte) *= sqrt(term); } } } @@ -433,20 +446,11 @@ void PhononThermalConductivity::calcFromRelaxons( } void PhononThermalConductivity::print() { - if (!mpi->mpiHead()) - return; - std::string units; - if (dimensionality == 1) { - units = "W m / K"; - } else if (dimensionality == 2) { - units = "W / K"; - } else { - units = "W / m / K"; - } + if (!mpi->mpiHead()) return; std::cout << "\n"; - std::cout << "Thermal Conductivity (" << units << ")\n"; + std::cout << "Thermal Conductivity (" << thCondUnits << ")\n"; for (int iCalc = 0; iCalc < statisticsSweep.getNumCalculations(); iCalc++) { auto calcStat = statisticsSweep.getCalcStatistics(iCalc); @@ -460,7 +464,7 @@ void PhononThermalConductivity::print() { std::cout << " " << std::scientific; for (int j = 0; j < dimensionality; j++) { std::cout << " " << std::setw(13) << std::right; - std::cout << tensordxd(iCalc, i, j) * thConductivityAuToSi; + std::cout << tensordxd(iCalc, i, j) * thCondConversion; } std::cout << "\n"; } @@ -469,17 +473,8 @@ void PhononThermalConductivity::print() { } void PhononThermalConductivity::outputToJSON(const std::string &outFileName) { - if (!mpi->mpiHead()) - return; - std::string units; - if (dimensionality == 1) { - units = "W m / K"; - } else if (dimensionality == 2) { - units = "W / K"; - } else { - units = "W /(m K)"; - } + if (!mpi->mpiHead()) return; std::vector temps; std::vector>> conductivities; @@ -495,7 +490,7 @@ void PhononThermalConductivity::outputToJSON(const std::string &outFileName) { for (int i = 0; i < dimensionality; i++) { std::vector cols; for (int j = 0; j < dimensionality; j++) { - cols.push_back(tensordxd(iCalc, i, j) * thConductivityAuToSi); + cols.push_back(tensordxd(iCalc, i, j) * thCondConversion); } rows.push_back(cols); } @@ -507,7 +502,7 @@ void PhononThermalConductivity::outputToJSON(const std::string &outFileName) { output["temperatures"] = temps; output["thermalConductivity"] = conductivities; output["temperatureUnit"] = "K"; - output["thermalConductivityUnit"] = units; + output["thermalConductivityUnit"] = thCondUnits; output["particleType"] = "phonon"; std::ofstream o(outFileName); o << std::setw(3) << output << std::endl; @@ -515,8 +510,8 @@ void PhononThermalConductivity::outputToJSON(const std::string &outFileName) { } void PhononThermalConductivity::print(const int &iter) { - if (!mpi->mpiHead()) - return; + + if (!mpi->mpiHead()) return; // get the time time_t currentTime; @@ -536,7 +531,7 @@ void PhononThermalConductivity::print(const int &iter) { std::cout.precision(5); for (int i = 0; i < dimensionality; i++) { std::cout << std::scientific; - std::cout << tensordxd(iCalc, i, i) * thConductivityAuToSi << " "; + std::cout << tensordxd(iCalc, i, i) * thCondConversion << " "; } std::cout << "\n"; } diff --git a/src/observable/phonon_thermal_cond.h b/src/observable/phonon_thermal_cond.h index 8f974e22..e4a23a1a 100644 --- a/src/observable/phonon_thermal_cond.h +++ b/src/observable/phonon_thermal_cond.h @@ -97,6 +97,10 @@ class PhononThermalConductivity : public Observable { protected: int whichType() override; BaseBandStructure &bandStructure; + // unit information for writing to files + std::string thCondUnits; + double thCondConversion; + }; #endif diff --git a/src/observable/phonon_viscosity.cpp b/src/observable/phonon_viscosity.cpp index bea4a093..c0a4e1c9 100644 --- a/src/observable/phonon_viscosity.cpp +++ b/src/observable/phonon_viscosity.cpp @@ -126,13 +126,13 @@ void PhononViscosity::calcFromRelaxons(Eigen::VectorXd &eigenvalues, // and stored ready to use here // calculate the first part of w^j_i,alpha - Eigen::Tensor tmpDriftEigvecs(3, 3, numStates); + Eigen::Tensor tmpDriftEigvecs(dimensionality, dimensionality, numStates); tmpDriftEigvecs.setZero(); for (int is : bandStructure.parallelStateIterator()) { auto isIdx = StateIndex(is); auto v = bandStructure.getGroupVelocity(isIdx); - for (int i : {0, 1, 2}) { - for (int j : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { + for (int j = 0; j < dimensionality; j++) { tmpDriftEigvecs(i, j, is) = phi(j, is) * v(i); } } @@ -140,10 +140,10 @@ void PhononViscosity::calcFromRelaxons(Eigen::VectorXd &eigenvalues, mpi->allReduceSum(&tmpDriftEigvecs); // now we're calculating w - Eigen::Tensor w(3, 3, numStates); + Eigen::Tensor w(dimensionality, dimensionality, numStates); w.setZero(); - for (auto i : {0, 1, 2}) { - for (auto j : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { + for (int j = 0; j < dimensionality; j++) { // drift eigenvectors * v -- only have phonon state indices // and cartesian directions Eigen::VectorXd x(numStates); diff --git a/src/observable/viscosity_io.cpp b/src/observable/viscosity_io.cpp index 8fffa8dd..d2868adf 100644 --- a/src/observable/viscosity_io.cpp +++ b/src/observable/viscosity_io.cpp @@ -7,10 +7,10 @@ // here alpha0 and alpha e are set through passing by reference void genericRelaxonEigenvectorsCheck(ParallelMatrix& eigenvectors, - int& numRelaxons, Particle& particle, - Eigen::VectorXd& theta0, - Eigen::VectorXd& theta_e, - int& alpha0, int& alpha_e) { + int& numRelaxons, Particle& particle, + Eigen::VectorXd& theta0, + Eigen::VectorXd& theta_e, + int& alpha0, int& alpha_e) { // calculate the overlaps with special eigenvectors Eigen::VectorXd prodTheta0(numRelaxons); prodTheta0.setZero(); @@ -19,6 +19,7 @@ void genericRelaxonEigenvectorsCheck(ParallelMatrix& eigenvectors, auto is = std::get<0>(tup); auto gamma = std::get<1>(tup); + //if(std::isnan(eigenvectors(is,gamma)) || std::isnan(theta0(is))) std::cout << is << " " << gamma << " " << eigenvectors(is,gamma) << " " << theta0(is) << std::endl; prodTheta0(gamma) += eigenvectors(is,gamma) * theta0(is); prodThetae(gamma) += eigenvectors(is,gamma) * theta_e(is); @@ -59,12 +60,12 @@ void genericRelaxonEigenvectorsCheck(ParallelMatrix& eigenvectors, // calculate special eigenvectors void genericCalcSpecialEigenvectors(BaseBandStructure& bandStructure, - StatisticsSweep& statisticsSweep, - double& spinFactor, - Eigen::VectorXd& theta0, - Eigen::VectorXd& theta_e, - Eigen::MatrixXd& phi, - double& C, Eigen::Vector3d& A) { + StatisticsSweep& statisticsSweep, + double& spinFactor, + Eigen::VectorXd& theta0, + Eigen::VectorXd& theta_e, + Eigen::MatrixXd& phi, + double& C, Eigen::Vector3d& A) { int dimensionality = bandStructure.getPoints().getCrystal().getDimensionality(); double volume = bandStructure.getPoints().getCrystal().getVolumeUnitCell(dimensionality); @@ -108,6 +109,7 @@ void genericCalcSpecialEigenvectors(BaseBandStructure& bandStructure, double U = 0; // specific heat + C = 0.; // calculate the special eigenvectors ---------------- for (int is : bandStructure.parallelStateIterator()) { @@ -123,6 +125,7 @@ void genericCalcSpecialEigenvectors(BaseBandStructure& bandStructure, theta_e(is) = sqrt(pop) * ds(is); U += pop; } + auto popCont = pop * (en - chemPot) * (en - chemPot); C += pop * (en - chemPot) * (en - chemPot); } mpi->allReduceSum(&theta0); @@ -131,13 +134,12 @@ void genericCalcSpecialEigenvectors(BaseBandStructure& bandStructure, mpi->allReduceSum(&U); // apply normalizations - C *= spinFactor / (volume * Npts * kBT * T); - theta0 *= 1./sqrt(kBT * T * volume * Npts * C); + C *= spinFactor / (volume * size_t(Npts) * kBT * T); + theta0 *= 1./sqrt(kBT * T * volume * size_t(Npts) * C); U *= spinFactor / (volume * Npts * kBT); if(particle.isPhonon()) U = 1.; // avoid making theta_e nan instead of zero theta_e *= 1./sqrt(kBT * U * Npts * volume); - - + // calculate A_i ---------------------------------------- // normalization coeff A ("phonon specific momentum") @@ -170,14 +172,14 @@ void genericCalcSpecialEigenvectors(BaseBandStructure& bandStructure, double pop = particle.getPopPopPm1(en, kBT, chemPot); // = n(n+1) auto q = bandStructure.getWavevector(isIdx); q = bandStructure.getPoints().bzToWs(q,Points::cartesianCoordinates); - for (auto i : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { phi(i, is) = q(i) * sqrt(pop) * ds(is); } } mpi->allReduceSum(&phi); // apply normalization to phi for(int is = 0; is < numStates; is++) { - for(int i : {0,1,2}) phi(i,is) *= 1./sqrt(kBT * volume * Npts * A(i)); + for (int i = 0; i < dimensionality; i++) phi(i,is) *= 1./sqrt(kBT * volume * Npts * A(i)); } /* // print phi overlap @@ -198,16 +200,26 @@ void printViscosity(std::string& viscosityName, Eigen::Tensor& viscos if (!mpi->mpiHead()) return; - std::string units; - if (dimensionality == 1) { units = "Pa s / m^2"; } // 3d - else if (dimensionality == 2) { units = "Pa s / m"; } // 2d - else { units = "Pa s"; } // 1d + // TODO Very important, do we have to multiply this by height / thickness? + std::string units, printHeader; + if (dimensionality == 1) { + units = "Pa s / m^2"; + printHeader = "i, eta[i,0]\n"; + } // 1d + else if (dimensionality == 2) { + units = "Pa s / m"; + printHeader = "i, j, eta[i,j,0], eta[i,j,1]\n"; + } // 2d + else { + units = "Pa s"; + printHeader = "i, j, k, eta[i,j,k,0], eta[i,j,k,1], eta[i,j,k,2]\n"; + } // 3d int numCalculations = statisticsSweep.getNumCalculations(); std::cout << "\n"; std::cout << viscosityName << " viscosity (" << units << ")\n"; - std::cout << "i, j, k, eta[i,j,k,0], eta[i,j,k,1], eta[i,j,k,2]\n"; + std::cout << printHeader; for (int iCalc = 0; iCalc < numCalculations; iCalc++) { @@ -222,7 +234,9 @@ void printViscosity(std::string& viscosityName, Eigen::Tensor& viscos for (int i = 0; i < dimensionality; i++) { for (int j = 0; j < dimensionality; j++) { for (int k = 0; k < dimensionality; k++) { - std::cout << i << " " << j << " " << k; + if(dimensionality == 1) std::cout << i; + if(dimensionality == 2) std::cout << i << " " << j; + if(dimensionality == 3) std::cout << i << " " << j << " " << k; for (int l = 0; l < dimensionality; l++) { std::cout << " " << std::setw(12) << std::right << viscosityTensor(iCalc, i, j, k, l) * viscosityAuToSi; @@ -243,7 +257,7 @@ void outputViscosityToJSON(const std::string& outFileName, const std::string& vi int numCalculations = statisticsSweep.getNumCalculations(); - std::string units; + std::string units; // TODO CHECK THIS! if (dimensionality == 1) { units = "Pa s / m^2"; } // 3d else if (dimensionality == 2) { units = "Pa s / m"; } // 2d else { units = "Pa s"; } // 1d @@ -317,18 +331,18 @@ void genericOutputRealSpaceToJSON(ScatteringMatrix& scatteringMatrix, auto calcStat = statisticsSweep.getCalcStatistics(0); // only one calc for relaxons double kBT = calcStat.temperature; - Eigen::MatrixXd Du(3,3); Du.setZero(); - Eigen::MatrixXd Wji0(3,3); Wji0.setZero(); + Eigen::MatrixXd Du(dimensionality,dimensionality); Du.setZero(); + Eigen::MatrixXd Wji0(dimensionality,dimensionality); Wji0.setZero(); // below used only for electrons - Eigen::MatrixXd Wjie(3,3); Wjie.setZero(); + Eigen::MatrixXd Wjie(dimensionality,dimensionality); Wjie.setZero(); // sum over the alpha and v states that this process owns for (auto tup : scatteringMatrix.getAllLocalStates()) { auto is1 = std::get<0>(tup); auto is2 = std::get<1>(tup); - for (auto i : {0, 1, 2}) { - for (auto j : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { + for (int j = 0; j < dimensionality; j++) { /* if(context.getUseUpperTriangle()) { if( i == j ) { Du(i,j) += phi(i,is1) * scatteringMatrix(is1,is2) * phi(j,is2); @@ -349,8 +363,8 @@ void genericOutputRealSpaceToJSON(ScatteringMatrix& scatteringMatrix, // discard acoustic phonon modes if (isPhonon && en < 0.001 / ryToCmm1) { continue; } auto v = bandStructure.getGroupVelocity(isIdx); - for (auto j : {0, 1, 2}) { - for (auto i : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { + for (int j = 0; j < dimensionality; j++) { // note: phi and theta here are elStates long, so we need to shift the state // index to account for the fact that we summed over the electronic part above // calculate qunatities for the real-space solve @@ -367,9 +381,9 @@ void genericOutputRealSpaceToJSON(ScatteringMatrix& scatteringMatrix, std::vector> vecDu; std::vector> vecWji0; std::vector> vecWjie; - for (auto i : {0, 1, 2}) { + for (int i = 0; i < dimensionality; i++) { std::vector temp1,temp2,temp3; - for (auto j : {0, 1, 2}) { + for (int j = 0; j < dimensionality; j++) { temp1.push_back(Du(i,j) / (energyRyToFs / twoPi)); temp2.push_back(Wji0(i,j) * velocityRyToSi); temp3.push_back(Wjie(i,j) * velocityRyToSi); @@ -425,7 +439,7 @@ void genericOutputRealSpaceToJSON(ScatteringMatrix& scatteringMatrix, output["specificHeatUnit"] = specificHeatUnits; output["particleType"] = particle.isPhonon() ? "phonon" : "electron"; std::vector Atemp; - for(int i = 0; i < 3; i++) { + for(int i = 0; i < dimensionality; i++) { Atemp.push_back(A(i) * Aconversion ); } output["Ai"] = Atemp; diff --git a/src/observable/wigner_phonon_thermal_cond.cpp b/src/observable/wigner_phonon_thermal_cond.cpp index 0f6bb610..c2badcdc 100644 --- a/src/observable/wigner_phonon_thermal_cond.cpp +++ b/src/observable/wigner_phonon_thermal_cond.cpp @@ -1,23 +1,34 @@ #include "wigner_phonon_thermal_cond.h" - #include - #include "constants.h" WignerPhononThermalConductivity::WignerPhononThermalConductivity( Context &context_, StatisticsSweep &statisticsSweep_, Crystal &crystal_, BaseBandStructure &bandStructure_, VectorBTE &relaxationTimes) - : PhononThermalConductivity(context_, statisticsSweep_, crystal_, - bandStructure_), + : PhononThermalConductivity(context_, statisticsSweep_, crystal_, bandStructure_), smaRelTimes(relaxationTimes) { - wignerCorrection = - Eigen::Tensor(numCalculations, dimensionality, dimensionality); + wignerCorrection = Eigen::Tensor(numCalculations, dimensionality, dimensionality); wignerCorrection.setZero(); auto particle = bandStructure.getParticle(); int dimensionality = crystal.getDimensionality(); + // set up units for writing to file + thCondUnits = "W /(m K)"; + if (dimensionality == 3) { + thCondConversion = thConductivityAuToSi; + } else if (dimensionality == 2) { + // multiply by the height of the cell / thickness of the cell to convert 3D -> 2D. + // Because the unit cell volume is already reduced for dimensionality, + // we only need to divide by thickness. + //double height = crystal.getDirectUnitCell()(2,2); + thCondConversion = thConductivityAuToSi * (1. / context.getThickness()); + } else { + Warning("1D conductivity should be manually adjusted for the cross-section of the material."); + thCondConversion = thConductivityAuToSi; + } + Eigen::VectorXd norm(numCalculations); for (int iCalc = 0; iCalc < numCalculations; iCalc++) { auto calcStat = statisticsSweep.getCalcStatistics(iCalc); @@ -153,20 +164,10 @@ void WignerPhononThermalConductivity::calcFromRelaxons( void WignerPhononThermalConductivity::print() { - if (!mpi->mpiHead()) - return; // debugging now - - std::string units; - if (dimensionality == 1) { - units = "W m / K"; - } else if (dimensionality == 2) { - units = "W / K"; - } else { - units = "W / m / K"; - } + if (!mpi->mpiHead()) return; std::cout << "\n"; - std::cout << "Wigner Thermal Conductivity (" << units << ")\n"; + std::cout << "Wigner Thermal Conductivity (" << thCondUnits << ")\n"; for (int iCalc = 0; iCalc < numCalculations; iCalc++) { auto calcStat = statisticsSweep.getCalcStatistics(iCalc); @@ -180,7 +181,7 @@ void WignerPhononThermalConductivity::print() { std::cout << " " << std::scientific; for (int j = 0; j < dimensionality; j++) { std::cout << " " << std::setw(13) << std::right; - std::cout << tensordxd(iCalc, i, j) * thConductivityAuToSi; + std::cout << tensordxd(iCalc, i, j) * thCondConversion; } std::cout << "\n"; } diff --git a/src/observable/wigner_phonon_thermal_cond.h b/src/observable/wigner_phonon_thermal_cond.h index 84d70d4f..233fac43 100644 --- a/src/observable/wigner_phonon_thermal_cond.h +++ b/src/observable/wigner_phonon_thermal_cond.h @@ -72,6 +72,9 @@ class WignerPhononThermalConductivity : public PhononThermalConductivity { protected: VectorBTE &smaRelTimes; Eigen::Tensor wignerCorrection; + // unit information for writing to files + std::string thCondUnits; + double thCondConversion; }; #endif diff --git a/src/parser/qe_input_parser.cpp b/src/parser/qe_input_parser.cpp index c6710dbf..92fcdad9 100644 --- a/src/parser/qe_input_parser.cpp +++ b/src/parser/qe_input_parser.cpp @@ -481,8 +481,8 @@ std::tuple QEParser::parsePhHarmonic(Context &context) { hasDielectric = true; } - // if there are the dielectric info, we can read dielectric matrix - // and the Born charges + // if there are the dielectric info, we can read dielectric matrix + // and the Born charges Eigen::Matrix3d dielectricMatrix = Eigen::Matrix3d::Zero(); Eigen::Tensor bornCharges(numAtoms, 3, 3); bornCharges.setZero(); @@ -656,7 +656,7 @@ QEParser::parseElHarmonicFourier(Context &context) { directUnitCell(1, 2) = std::stod(lineSplit[1]); directUnitCell(2, 2) = std::stod(lineSplit[2]); - // Now we parse the electronic structure + // Now we parse the electronic structure ============================ pugi::xml_node bandStructureXML = output.child("band_structure"); bool isLSDA = bandStructureXML.child("lsda").text().as_bool(); @@ -666,7 +666,7 @@ QEParser::parseElHarmonicFourier(Context &context) { // note: nelec is written as double in the XML file! int numElectrons = int(bandStructureXML.child("nelec").text().as_double()); - // get fermi energy + // get fermi energy or HOMO --------------------- double homo; // it's an insulator if (bandStructureXML.child("highestOccupiedLevel")) { @@ -678,6 +678,8 @@ QEParser::parseElHarmonicFourier(Context &context) { "nor fermi_energy tags appear in XML file."); } homo *= 2.;// conversion from Hartree to Rydberg + + // parse the points grid ------------------------------------- int numIrreduciblePoints = bandStructureXML.child("nks").text().as_int(); pugi::xml_node startingKPoints = bandStructureXML.child("starting_k_points"); @@ -690,16 +692,15 @@ QEParser::parseElHarmonicFourier(Context &context) { "doesn't make sense for Phoebe -- likely you forgot to perform NSCF on\n" "the full k-mesh first."); } - // Error("Grid found in QE:XML, should have used full kPoints grid"); + Error("Grid found in QE:XML, should have used full kPoints grid"); } // Initialize the crystal class - Crystal crystal(context, directUnitCell, atomicPositions, atomicSpecies, speciesNames, speciesMasses); crystal.print(); - // initialize reciprocal lattice cell + // initialize reciprocal lattice cell -------------------------- // I need this to convert kPoints from cartesian to crystal coordinates pugi::xml_node basisSet = output.child("basis_set"); @@ -718,7 +719,7 @@ QEParser::parseElHarmonicFourier(Context &context) { bVectors(1, 2) = std::stod(lineSplit[1]); bVectors(2, 2) = std::stod(lineSplit[2]); - // parse k-points and energies + // parse k-points and energies ---------------------------- Eigen::Matrix irrPoints(3, numIrreduciblePoints); Eigen::VectorXd irrWeights(numIrreduciblePoints); diff --git a/src/statistics_sweep.cpp b/src/statistics_sweep.cpp index 5e97e54b..c8e92258 100644 --- a/src/statistics_sweep.cpp +++ b/src/statistics_sweep.cpp @@ -277,8 +277,8 @@ StatisticsSweep::findChemicalPotentialFromDoping(const double &doping, // numElectrons is the number of electrons in the unit cell before doping // doping > 0 means p-doping (fermi level in the valence band) - numElectronsDoped = - occupiedStates - doping * volume * pow(distanceBohrToCm, 3) / spinFactor; + numElectronsDoped = occupiedStates - doping * volume * pow(distanceBohrToCm, 3) / spinFactor; + // bisection method: I need to find the root of N - \int fermi dirac = 0 // initial guess double chemicalPotential = fermiLevel; @@ -295,18 +295,18 @@ StatisticsSweep::findChemicalPotentialFromDoping(const double &doping, "were not excluded using exclude_bands in Wannier90.\nSee a note about this in the elphWannier tutorial."); } if (numElectronsDoped < 0.) { - Error("The number of occupied states is negative"); + Error("The number of occupied states is negative!"); } // if we are looking for the fermi level at T=0 and n=0, we have a corner case // when we have completely empty bands or completely full bands. if (doping == 0. && temperature == 0.) {// case of computing fermi level if (numElectronsDoped == 0.) { - fermiLevel = *min_element(energies.begin(), energies.end()); + if(energies.size() > 0) fermiLevel = *min_element(energies.begin(), energies.end()); chemicalPotential = fermiLevel; return chemicalPotential; } else if (numElectronsDoped == float(numBands)) { - fermiLevel = *max_element(energies.begin(), energies.end()); + if(energies.size() > 0) fermiLevel = *max_element(energies.begin(), energies.end()); chemicalPotential = fermiLevel; return chemicalPotential; } @@ -334,12 +334,12 @@ StatisticsSweep::findChemicalPotentialFromDoping(const double &doping, // check if starting values are bad if (sgn(aY) == sgn(bY)) { - Error("I should revisit the boundary limits for bisection method"); + Error("Something is wrong with the boundary limits for mu determination bisection method."); } for (int iter = 0; iter < maxIter; iter++) { if (mpi->mpiHead() && iter == maxIter - 1) { - Error("Max iteration reached in finding mu"); + Error("Max iteration reached without finding mu."); } // x value is midpoint of prior values double cX = (aX + bX) / 2.; @@ -415,8 +415,7 @@ void StatisticsSweep::printInfo() { std::cout << "Statistical parameters for the calculation\n"; if (particle.isElectron()) { - std::cout << "Fermi level: " << fermiLevel * energyRyToEv << " (eV)" - << std::endl; + std::cout << "Fermi level: " << fermiLevel * energyRyToEv << " (eV)" << std::endl; std::cout << "Index, temperature, chemical potential, doping concentration\n"; } else { std::cout << "Index, temperature\n"; diff --git a/src/utilities.cpp b/src/utilities.cpp index 6fc05eaa..3dd4c511 100644 --- a/src/utilities.cpp +++ b/src/utilities.cpp @@ -79,9 +79,11 @@ std::tuple memoryUsage() { double findMaxRelativeDifference(const Eigen::Tensor &x, const Eigen::Tensor &xRef) { + if (!(x.dimensions() == xRef.dimensions())) { - Error("Can't compare inconsistent tensors"); + Error("Developer error: Can't compare inconsistent tensors"); } + Eigen::Tensor diff = ((x - xRef) / xRef).abs(); double maxDiff = 1.e20; for (int i = 0; i < diff.dimension(0); i++) { @@ -94,7 +96,7 @@ double findMaxRelativeDifference(const Eigen::Tensor &x, } } if (maxDiff == 1.e20) { - Error("wrong bounds in findMaxRelativeDifference"); + Error("Developer error: Wrong bounds in findMaxRelativeDifference."); } return maxDiff; }