1 /*------------------------------------------------------------------- 2 Copyright 2011 Ravishankar Sundararaman 3 4 This file is part of JDFTx. 5 6 JDFTx is free software: you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation, either version 3 of the License, or 9 (at your option) any later version. 10 11 JDFTx is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>. 18 -------------------------------------------------------------------*/ 19 20 #include <commands/command.h> 21 #include <electronic/Everything.h> 22 #include <electronic/ColumnBundle.h> 23 #include <electronic/Polarizability.h> 24 #include <electronic/ElectronScattering.h> 25 #include <electronic/Vibrations.h> 26 #include <electronic/Dump_internal.h> 27 #include <electronic/DumpBGW_internal.h> 28 #include <core/Units.h> 29 30 struct CommandDumpOnly : public Command 31 { CommandDumpOnlyCommandDumpOnly32 CommandDumpOnly() : Command("dump-only", "jdftx/Output") 33 { 34 comments = 35 "Bypass all minimization, perform a single energy evaluation at\n" 36 "the initial state, and process dump commands at the end state.\n" 37 "Useful for post-processing a converged calculation."; 38 39 forbid("fix-electron-potential"); 40 forbid("fix-electron-density"); 41 } 42 processCommandDumpOnly43 void process(ParamList& pl, Everything& e) 44 { e.cntrl.dumpOnly = true; 45 } 46 printStatusCommandDumpOnly47 void printStatus(Everything& e, int iRep) 48 { 49 } 50 } 51 commandDumpOnly; 52 53 //Variables allowed at dump frequency Init 54 EnumStringMap<DumpVariable> varInitMap 55 ( DumpNone, "None", 56 DumpIonicPositions, "IonicPositions", 57 DumpLattice, "Lattice", 58 DumpCoreDensity, "CoreDensity", 59 DumpVlocps, "Vlocps", 60 DumpSymmetries, "Symmetries", 61 DumpKpoints, "Kpoints", 62 DumpGvectors, "Gvectors" 63 ); 64 65 EnumStringMap<DumpFrequency> freqMap 66 ( DumpFreq_End, "End", 67 DumpFreq_Init, "Init", 68 DumpFreq_Ionic, "Ionic", 69 DumpFreq_Gummel, "Gummel", 70 DumpFreq_Fluid, "Fluid", 71 DumpFreq_Electronic, "Electronic" 72 ); 73 EnumStringMap<DumpFrequency> freqDescMap 74 ( DumpFreq_End, "Dump specified vars at the end of the calculation", 75 DumpFreq_Init, "Dump specified vars from " + varInitMap.optionList() + " after initialization (even in dry run)", 76 DumpFreq_Ionic, "Dump specified vars every (few) ionic / lattice step(s)", 77 DumpFreq_Gummel, "Dump specified vars every (few) fluid+electron minimize of the gummel loop", 78 DumpFreq_Fluid, "Dump specified vars every (few) fluid step(s)", 79 DumpFreq_Electronic, "Dump specified vars every (few) electronic step(s)" 80 ); 81 82 83 EnumStringMap<DumpVariable> varMap 84 ( DumpNone, "None", 85 DumpState, "State", 86 DumpIonicPositions, "IonicPositions", 87 DumpForces, "Forces", 88 DumpLattice, "Lattice", 89 DumpIonicDensity, "IonicDensity", 90 DumpElecDensity, "ElecDensity", 91 DumpElecDensityAccum, "ElecDensityAccum", 92 DumpCoreDensity, "CoreDensity", 93 DumpKEdensity, "KEdensity", 94 DumpFluidDensity, "FluidDensity", 95 DumpDvac, "Dvac", 96 DumpDfluid, "Dfluid", 97 DumpDtot, "Dtot", 98 DumpVcavity, "Vcavity", 99 DumpVfluidTot, "VfluidTot", 100 DumpVlocps, "Vlocps", 101 DumpVscloc, "Vscloc", 102 DumpBandEigs, "BandEigs", 103 DumpBandProjections, "BandProjections", 104 DumpEigStats, "EigStats", 105 DumpFillings, "Fillings", 106 DumpRhoAtom, "RhoAtom", 107 DumpBandUnfold, "BandUnfold", 108 DumpEcomponents, "Ecomponents", 109 DumpExcCompare, "ExcCompare", 110 DumpBoundCharge, "BoundCharge", 111 DumpSolvationRadii, "SolvationRadii", 112 DumpQMC, "QMC", 113 DumpOcean, "Ocean", 114 DumpBGW, "BGW", 115 DumpRealSpaceWfns, "RealSpaceWfns", 116 DumpFluidDebug, "FluidDebug", 117 DumpSlabEpsilon, "SlabEpsilon", 118 DumpBulkEpsilon, "BulkEpsilon", 119 DumpChargedDefect, "ChargedDefect", 120 DumpDOS, "DOS", 121 DumpSIC, "SelfInteractionCorrection", 122 DumpDipole, "Dipole", 123 DumpStress, "Stress", 124 DumpExcitations, "Excitations", 125 DumpFCI, "FCI", 126 DumpSpin, "Spin", 127 DumpMomenta, "Momenta", 128 DumpVelocities, "Velocities", 129 DumpFermiVelocity, "FermiVelocity", 130 DumpSymmetries, "Symmetries", 131 DumpKpoints, "Kpoints", 132 DumpGvectors, "Gvectors", 133 DumpOrbitalDep, "OrbitalDep", 134 DumpXCanalysis, "XCanalysis", 135 DumpEresolvedDensity, "EresolvedDensity", 136 DumpFermiDensity, "FermiDensity" 137 ); 138 EnumStringMap<DumpVariable> varDescMap 139 ( DumpNone, "Dump nothing", 140 DumpState, "All variables needed to restart calculation: wavefunction and fluid state/fillings if any", 141 DumpIonicPositions, "Ionic positions in the same format (and coordinate system) as the input file", 142 DumpForces, "Forces on the ions in the coordinate system selected by command forces-output-coords", 143 DumpLattice, "Lattice vectors in the same format as the input file", 144 DumpIonicDensity, "Nuclear charge density (with gaussians)", 145 DumpElecDensity, "Electronic densities (n or nup,ndn)", 146 DumpElecDensityAccum, "Electronic densities (n or nup,ndn) accumulated over MD trajectory", 147 DumpCoreDensity, "Total core electron density (from partial core corrections)", 148 DumpKEdensity, "Kinetic energy density of the valence electrons", 149 DumpFluidDensity, "Fluid densities (NO,NH,nWater for explicit fluids, cavity function for PCMs)", 150 DumpDvac, "Electrostatic potential due to explicit system alone", 151 DumpDfluid, "Electrostatic potential due to fluid alone", 152 DumpDtot, "Total electrostatic potential", 153 DumpVcavity, "Fluid cavitation potential on the electron density that determines the cavity", 154 DumpVfluidTot, "Total contribution of fluid to the electron potential", 155 DumpVlocps, "Local part of pseudopotentials", 156 DumpVscloc, "Self-consistent potential", 157 DumpBandEigs, "Band Eigenvalues", 158 DumpBandProjections,"Projections of each band state against each atomic orbital", 159 DumpEigStats, "Band eigenvalue statistics: HOMO, LUMO, min, max and Fermi level", 160 DumpFillings, "Fillings", 161 DumpRhoAtom, "Atomic-orbital projected density matrices (only for species with +U enabled)", 162 DumpBandUnfold, "Unfold band structure from supercell to unit cell (see command band-unfold)", 163 DumpEcomponents, "Components of the energy", 164 DumpBoundCharge, "Bound charge in the fluid", 165 DumpSolvationRadii, "Effective solvation radii based on fluid bound charge distribution", 166 DumpQMC, "Blip'd orbitals and potential for CASINO \\cite Katie-QMC", 167 DumpOcean, "Wave functions for Ocean code", 168 DumpBGW, "G-space wavefunctions, density and potential for Berkeley GW (requires HDF5 support)", 169 DumpRealSpaceWfns, "Real-space wavefunctions (one column per file)", 170 DumpExcCompare, "Energies for other exchange-correlation functionals (see command elec-ex-corr-compare)", 171 DumpFluidDebug, "Fluid specific debug output if any ", 172 DumpSlabEpsilon, "Local dielectric function of a slab (see command slab-epsilon) ", 173 DumpBulkEpsilon, "Dielectric constant of a periodic solid (see command bulk-epsilon) ", 174 DumpChargedDefect, "Calculate energy correction for charged defect (see command charged-defect) ", 175 DumpDOS, "Density of States (see command density-of-states)", 176 DumpSIC, "Calculates Perdew-Zunger self-interaction corrected Kohn-Sham eigenvalues", 177 DumpDipole, "Dipole moment of explicit charges (ionic and electronic)", 178 DumpStress, "Dumps dE/dR_ij where R_ij is the i'th component of the j'th lattice vector", 179 DumpExcitations, "Dumps dipole moments and transition strength (electric-dipole) of excitations", 180 DumpFCI, "Output Coulomb matrix elements in FCIDUMP format", 181 DumpSpin, "Spin matrix elements from non-collinear calculations in a binary file (indices outer to inner: state, cartesian direction, band1, band2)", 182 DumpMomenta, "Momentum matrix elements in a binary file (indices outer to inner: state, cartesian direction, band1, band2)", 183 DumpVelocities, "Diagonal momentum/velocity matrix elements in a binary file (indices outer to inner: state, band, cartesian direction)", 184 DumpFermiVelocity, "Fermi velocity, density of states at Fermi level and related quantities", 185 DumpSymmetries, "List of symmetry matrices (in covariant lattice coordinates)", 186 DumpKpoints, "List of reduced k-points in calculation, and mapping to the unreduced k-point mesh", 187 DumpGvectors, "List of G vectors in reciprocal lattice basis, for each k-point", 188 DumpOrbitalDep, "Custom output from orbital-dependent functionals (eg. quasi-particle energies, discontinuity potential)", 189 DumpXCanalysis, "Debug VW KE density, single-particle-ness and spin-polarzied Hartree potential", 190 DumpEresolvedDensity, "Electron density from bands within specified energy ranges", 191 DumpFermiDensity, "Electron density from fermi-derivative at specified energy" 192 ); 193 194 struct CommandDump : public Command 195 { CommandDumpCommandDump196 CommandDump() : Command("dump", "jdftx/Output") 197 { 198 format = "<freq> <var> <var> ..."; 199 comments = 200 "<freq> is one of:" 201 + addDescriptions(freqMap.optionList(), linkDescription(freqMap, freqDescMap)) 202 + "\n\nand each <var> is one of:" 203 + addDescriptions(varMap.optionList(), linkDescription(varMap, varDescMap)) 204 + "\n\nList of dumped variables from multiple instances will be accumulated for each <freq>." 205 "\nUse command dump-interval to dump at regular intervals instead of every iteration."; 206 allowMultiple = true; 207 hasDefault = true; 208 } 209 processCommandDump210 void process(ParamList& pl, Everything& e) 211 { DumpFrequency freq; 212 pl.get(freq, DumpFreq_Delim, freqMap, "freq"); 213 //Handle the default call 214 if(freq==DumpFreq_Delim) 215 { e.dump.insert(std::make_pair(DumpFreq_End,DumpState)); 216 return; 217 } 218 //For any real dump frequency: 219 while(true) 220 { DumpVariable var; 221 pl.get(var, DumpDelim, freq==DumpFreq_Init ? varInitMap : varMap, "var"); 222 if(var==DumpDelim) break; //will happen at end of command line 223 e.dump.insert(std::make_pair(freq,var)); 224 //Check for unsupported features: 225 #ifndef HDF5_ENABLED 226 if(var==DumpBGW) throw string("BerkeleyGW interface requires HDF5 support (CMake option EnableHDF5)\n"); 227 #endif 228 } 229 } 230 printStatusCommandDump231 void printStatus(Everything& e, int iRep) 232 { //Coealesce dump outputs for each frequency into a single command 233 std::multimap<DumpFrequency,DumpVariable> dumpMap(e.dump.begin(), e.dump.end()); 234 typedef std::multimap<DumpFrequency,DumpVariable>::iterator Iter; 235 Iter i = dumpMap.begin(); 236 int iDump=0; 237 while(i!=dumpMap.end()) 238 { if(iDump==iRep) logPrintf("%s", freqMap.getString(i->first)); 239 std::pair<Iter,Iter> range = dumpMap.equal_range(i->first); 240 for(i=range.first; i!=range.second; i++) 241 if(iDump==iRep) logPrintf(" %s", varMap.getString(i->second)); 242 iDump++; 243 } 244 } 245 } 246 commandDump; 247 248 249 struct CommandDumpInterval : public Command 250 { CommandDumpIntervalCommandDumpInterval251 CommandDumpInterval() : Command("dump-interval", "jdftx/Output") 252 { 253 format = "<freq> <interval>"; 254 comments = 255 "Dump every <interval> iterations of type <freq>=Ionic|Electronic|Fluid|Gummel.\n" 256 "Without this command, the behavior defaults to <interval>=1 for each <freq>.\n" 257 "(See command dump for more details)"; 258 allowMultiple = true; 259 } 260 processCommandDumpInterval261 void process(ParamList& pl, Everything& e) 262 { //get the frequency: 263 DumpFrequency freq; 264 pl.get(freq, DumpFreq_Delim, freqMap, "freq", true); 265 if(freq==DumpFreq_End || freq==DumpFreq_Init) 266 throw string("<freq> must be one of Ionic|Electronic|Fluid|Gummel"); 267 if(e.dump.interval.find(freq) != e.dump.interval.end()) 268 throw string("dump-interval has been specified multiple times for <freq>=") + freqMap.getString(freq); 269 //get the interval: 270 int interval; 271 pl.get(interval, 1, "interval", true); 272 if(interval<1) 273 throw string("<interval> must be a positive integer"); 274 //Set the interval 275 e.dump.interval[freq] = interval; 276 } 277 printStatusCommandDumpInterval278 void printStatus(Everything& e, int iRep) 279 { std::map<DumpFrequency,int>::const_iterator iter = e.dump.interval.begin(); 280 for(int i=0; i<iRep; i++) iter++; //access the iRep'th entry of interval 281 logPrintf("%s %d", freqMap.getString(iter->first), iter->second); 282 } 283 } 284 commandDumpInterval; 285 286 287 struct CommandDumpName : public Command 288 { CommandDumpNameCommandDumpName289 CommandDumpName() : Command("dump-name", "jdftx/Output") 290 { 291 format = "<format> [<freq1> <format1>] [<freq2> <format2>] ..."; 292 comments = 293 "Control the filename pattern for dump output, where <format> is an\n" 294 "arbitrary format string that will be substituted according to:\n" 295 "+ $VAR -> name of the variable being dumped (this must be present)\n" 296 "+ $ITER -> iteration number of relevant dump frequency\n" 297 "+ $INPUT -> base name of input file, or 'stdin'\n" 298 "+ $STAMP -> time-stamp at the start of dump\n" 299 "\n" 300 "Optionally, a different <format> could be specified for some dump frequencies."; 301 hasDefault = true; 302 } 303 processCommandDumpName304 void process(ParamList& pl, Everything& e) 305 { pl.get(e.dump.format, string("$INPUT.$VAR"), "format"); 306 if(e.dump.format.find("$VAR")==string::npos) 307 throw "<format> = " + e.dump.format + " doesn't contain the pattern $VAR"; 308 //Check for additional specifictaions: 309 while(true) 310 { DumpFrequency freq; pl.get(freq, DumpFreq_Delim, freqMap, "<freqN>"); 311 if(freq==DumpFreq_Delim) break; //no more freq's 312 string format; pl.get(format, string(), "<formatN>", true); 313 if(format.find("$VAR")==string::npos) 314 throw "<format> = " + format + " doesn't contain the pattern $VAR"; 315 e.dump.formatFreq[freq] = format; 316 } 317 } 318 printStatusCommandDumpName319 void printStatus(Everything& e, int iRep) 320 { logPrintf("%s", e.dump.format.c_str()); 321 for(auto entry: e.dump.formatFreq) 322 logPrintf(" \\\n\t%s %s", freqMap.getString(entry.first), entry.second.c_str()); 323 } 324 } 325 commandDumpName; 326 327 328 EnumStringMap<Polarizability::EigenBasis> polarizabilityMap 329 ( Polarizability::NonInteracting, "NonInteracting", 330 Polarizability::External, "External", 331 Polarizability::Total, "Total" 332 ); 333 334 struct CommandPolarizability : public Command 335 { CommandPolarizabilityCommandPolarizability336 CommandPolarizability() : Command("polarizability", "jdftx/Output") 337 { 338 format = "<eigenBasis>=" + polarizabilityMap.optionList() + " [<Ecut>=0] [<nEigs>=0]"; 339 comments = "Output polarizability matrix in specified eigenBasis"; 340 341 forbid("electron-scattering"); //both are major operations that are given permission to destroy Everything if necessary 342 } 343 processCommandPolarizability344 void process(ParamList& pl, Everything& e) 345 { e.dump.polarizability = std::make_shared<Polarizability>(); 346 pl.get(e.dump.polarizability->eigenBasis, Polarizability::NonInteracting, polarizabilityMap, "eigenBasis"); 347 pl.get(e.dump.polarizability->Ecut, 0., "Ecut"); 348 pl.get(e.dump.polarizability->nEigs, 0, "nEigs"); 349 e.dump.insert(std::make_pair(DumpFreq_End, DumpPolarizability)); 350 } 351 printStatusCommandPolarizability352 void printStatus(Everything& e, int iRep) 353 { logPrintf("%s %lg %d", polarizabilityMap.getString(e.dump.polarizability->eigenBasis), 354 e.dump.polarizability->Ecut, e.dump.polarizability->nEigs); 355 } 356 } 357 commandPolarizability; 358 359 360 enum ElectronScatteringMember 361 { ESM_eta, 362 ESM_Ecut, 363 ESM_fCut, 364 ESM_omegaMax, 365 ESM_RPA, 366 ESM_slabResponse, 367 ESM_EcutTransverse, 368 ESM_delim 369 }; 370 EnumStringMap<ElectronScatteringMember> esmMap 371 ( ESM_eta, "eta", 372 ESM_Ecut, "Ecut", 373 ESM_fCut, "fCut", 374 ESM_omegaMax, "omegaMax", 375 ESM_RPA, "RPA", 376 ESM_slabResponse, "slabResponse", 377 ESM_EcutTransverse, "EcutTransverse" 378 ); 379 380 struct CommandElectronScattering : public Command 381 { CommandElectronScatteringCommandElectronScattering382 CommandElectronScattering() : Command("electron-scattering", "jdftx/Output") 383 { 384 format = "<key1> <value1> ..."; 385 comments = "Calculate electron-electron scattering rates (expensive!)\n" 386 "and output contribution to imaginary part of electron self-energy\n" 387 "(calculated effectively using full-frequency G0W0).\n" 388 "\n" 389 "The following key-value pairs can appear in any order:\n" 390 "\n+ eta <eta>\n\n" 391 " <eta> in Eh specifies frequency grid resolution (required)\n" 392 "\n+ Ecut <Ecut>\n\n" 393 " <Ecut> in Eh specifies energy cut-off for dielectric matrices.\n" 394 " (If zero, the wavefunction cutoff from elec-cutoff is used instead.)\n" 395 "\n+ fCut <fCut>\n\n" 396 " <fCut> specifies threshold for considering states fully occupied or\n" 397 " unoccupied in optimizing sums over states (default: 1e-6)\n" 398 "\n+ omegaMax <omegaMax>\n\n" 399 " <omegaMax> in Eh is the maximum energy transfer to account for\n" 400 " and hence the maximum frequency in dielectric function frequency grid.\n" 401 " (if zero, autodetermine from available eigenvalues)\n" 402 "\n+ RPA yes|no\n\n" 403 " If yes, use RPA response that ignores XC contribution. (default: no).\n" 404 "\n+ slabResponse yes|no\n\n" 405 " Whether to output slab-normal-direction susceptibility instead.\n" 406 " This needs slab geometry in coulomb-interaction, and will bypass the\n" 407 " actual electron-electron scattering calculation and output.\n" 408 "\n+ EcutTransverse <EcutTransverse>\n\n" 409 " <EcutTransverse> in Eh specifies energy cut-off for dielectric matrix in.\n" 410 " directions trasverse to the slab normal; only valid when slabResponse = yes.\n" 411 " (If zero, use the same value as Ecut above.)"; 412 413 require("coulomb-interaction"); 414 forbid("polarizability"); //both are major operations that are given permission to destroy Everything if necessary 415 } 416 processCommandElectronScattering417 void process(ParamList& pl, Everything& e) 418 { e.dump.electronScattering = std::make_shared<ElectronScattering>(); 419 e.dump.insert(std::make_pair(DumpFreq_End, DumpElectronScattering)); 420 ElectronScattering& es = *(e.dump.electronScattering); 421 while(true) 422 { ElectronScatteringMember key; 423 pl.get(key, ESM_delim, esmMap, "key"); 424 if(key == ESM_delim) break; //end of input 425 switch(key) 426 { case ESM_eta: pl.get(es.eta, 0., "eta", true); break; 427 case ESM_Ecut: pl.get(es.Ecut, 0., "Ecut", true); break; 428 case ESM_fCut: pl.get(es.fCut, 0., "fCut", true); break; 429 case ESM_omegaMax: pl.get(es.omegaMax, 0., "omegaMax", true); break; 430 case ESM_RPA: pl.get(es.RPA, false, boolMap, "RPA", true); break; 431 case ESM_slabResponse: pl.get(es.slabResponse, false, boolMap, "slabResponse", true); break; 432 case ESM_EcutTransverse: pl.get(es.EcutTransverse, 0., "EcutTransverse", true); break; 433 case ESM_delim: break; //never encountered; to suppress compiler warning 434 } 435 } 436 if(es.slabResponse) 437 { if(e.coulombParams.geometry != CoulombParams::Slab) 438 throw string("slabResponse = yes requires slab geometry in coulomb-interaction"); 439 } 440 else 441 { if(es.EcutTransverse) throw string("Cannot specify EcutTransverse when slabResponse = no"); 442 } 443 if(es.eta <= 0.) throw string("Must specify frequency grid resolution eta > 0."); 444 } 445 printStatusCommandElectronScattering446 void printStatus(Everything& e, int iRep) 447 { const ElectronScattering& es = *(e.dump.electronScattering); 448 logPrintf(" \\\n\teta %lg", es.eta); 449 logPrintf(" \\\n\tEcut %lg", es.Ecut); 450 logPrintf(" \\\n\tfCut %lg", es.fCut); 451 logPrintf(" \\\n\tomegaMax %lg", es.omegaMax); 452 logPrintf(" \\\n\tRPA %s", boolMap.getString(es.RPA)); 453 logPrintf(" \\\n\tslabResponse %s", boolMap.getString(es.slabResponse)); 454 if(es.slabResponse) logPrintf(" \\\n\tEcutTransverse %lg", es.EcutTransverse); 455 } 456 } 457 commandElectronScattering; 458 459 460 struct CommandPolarizabilityKdiff : public Command 461 { CommandPolarizabilityKdiffCommandPolarizabilityKdiff462 CommandPolarizabilityKdiff() : Command("polarizability-kdiff", "jdftx/Output") 463 { 464 format = "<dk0> <dk1> <dk2> [<dkFilenamePattern>]"; 465 comments = "Select k-point difference (in reciprocal lattice coords) for polarizability output.\n" 466 "\n" 467 "<dkFilenamePattern> may be specified to read offset band structure calcualations when <dk>\n" 468 "does not belong to the k-point mesh. This string should be a filename pattern containing\n" 469 "$VAR (to be replaced by eigenvals and wfns) and $q (to be replaced by state index)"; 470 471 require("polarizability"); 472 } 473 processCommandPolarizabilityKdiff474 void process(ParamList& pl, Everything& e) 475 { pl.get(e.dump.polarizability->dk[0], 0., "dk0", true); 476 pl.get(e.dump.polarizability->dk[1], 0., "dk1", true); 477 pl.get(e.dump.polarizability->dk[2], 0., "dk2", true); 478 //Optional filename for offset states: 479 string& dkFilenamePattern = e.dump.polarizability->dkFilenamePattern; 480 pl.get(dkFilenamePattern, string(), "dkFilenamePattern"); 481 if(dkFilenamePattern.length()) 482 { if(dkFilenamePattern.find("$VAR")==string::npos) throw "<dkFilenamePattern> = " + dkFilenamePattern + " doesn't contain the pattern $VAR"; 483 if(dkFilenamePattern.find("$q")==string::npos) throw "<dkFilenamePattern> = " + dkFilenamePattern + " doesn't contain the pattern $q"; 484 } 485 } 486 printStatusCommandPolarizabilityKdiff487 void printStatus(Everything& e, int iRep) 488 { for(int i=0; i<3; i++) logPrintf("%lg ", e.dump.polarizability->dk[i]); 489 logPrintf("%s", e.dump.polarizability->dkFilenamePattern.c_str()); 490 } 491 } 492 commandPolarizabilityKdiff; 493 494 495 struct CommandDumpEresolvedDensity : public Command 496 { CommandDumpEresolvedDensityCommandDumpEresolvedDensity497 CommandDumpEresolvedDensity() : Command("dump-Eresolved-density", "jdftx/Output") 498 { 499 format = "<Emin> <Emax>"; 500 comments = 501 "Output electron density from bands within a specified energy range\n" 502 "[Emin,Emax] (in Eh).\n" 503 "\n" 504 "When issued multiple times, the outputs will be\n" 505 "numbered sequenetially EresolvedDensity.0 etc.\n" 506 "This automatically invokes dump at End; dumping at\n" 507 "other frequencies may be requested using the dump command."; 508 509 allowMultiple = true; 510 } 511 processCommandDumpEresolvedDensity512 void process(ParamList& pl, Everything& e) 513 { double Emin, Emax; 514 pl.get(Emin, 0., "Emin", true); 515 pl.get(Emax, 0., "Emax", true); 516 if(Emin >= Emax) throw string("Emin must be < Emax"); 517 e.dump.densityErange.push_back(std::make_pair(Emin, Emax)); 518 e.dump.insert(std::make_pair(DumpFreq_End, DumpEresolvedDensity)); 519 } 520 printStatusCommandDumpEresolvedDensity521 void printStatus(Everything& e, int iRep) 522 { const auto& Erange = e.dump.densityErange[iRep]; 523 logPrintf("%lg %lg", Erange.first, Erange.second); 524 } 525 } 526 commandDumpEresolvedDensity; 527 528 struct CommandDumpFermiDensity : public Command 529 { CommandDumpFermiDensityCommandDumpFermiDensity530 CommandDumpFermiDensity() : Command("dump-fermi-density", "jdftx/Output") 531 { 532 format = "[<muLevel>]"; 533 comments = 534 "Output electron density calculated from derivative of smearing\n" 535 "function evaluated at desired chemical potential <muLevel>.\n" 536 "If unspecified, calculated chemical potential will be used.\n" 537 "\n" 538 "When issued multiple times, the outputs will be\n" 539 "numbered sequenetially FermiDensity.0 etc.\n" 540 "This automatically invokes dump at End; dumping at\n" 541 "other frequencies may be requested using the dump command."; 542 543 allowMultiple = true; 544 require("elec-smearing"); //require smearing 545 } 546 processCommandDumpFermiDensity547 void process(ParamList& pl, Everything& e) 548 { double muLevel; 549 pl.get(muLevel, std::numeric_limits<double>::quiet_NaN(), "muLevel", false); 550 e.dump.fermiDensityLevels.push_back(muLevel); 551 e.dump.insert(std::make_pair(DumpFreq_End, DumpFermiDensity)); 552 } 553 printStatusCommandDumpFermiDensity554 void printStatus(Everything& e, int iRep) 555 { const auto& muLevel = e.dump.fermiDensityLevels[iRep]; 556 logPrintf("%lg", muLevel); 557 } 558 } 559 commandDumpFermiDensity; 560 561 562 //---------- Vibrations ------------------------ 563 564 //An enum corresponding to members of class Vibrations 565 enum VibrationsMember 566 { VM_dr, 567 VM_centralDiff, 568 VM_useConstraints, 569 VM_translationSym, 570 VM_rotationSym, 571 VM_omegaMin, 572 VM_T, 573 VM_omegaResolution, 574 VM_Delim 575 }; 576 577 EnumStringMap<VibrationsMember> vibMap 578 ( VM_dr, "dr", 579 VM_centralDiff, "centralDiff", 580 VM_useConstraints, "useConstraints", 581 VM_translationSym, "translationSym", 582 VM_rotationSym, "rotationSym", 583 VM_omegaMin, "omegaMin", 584 VM_T, "T", 585 VM_omegaResolution, "omegaResolution" 586 ); 587 588 struct CommandVibrations : public Command 589 { CommandVibrationsCommandVibrations590 CommandVibrations() : Command("vibrations", "jdftx/Output") 591 { 592 format = "<key1> <args1> ..."; 593 comments = 594 "Calculate vibrational modes of the system using a finite difference method.\n" 595 "Note that this command should typically be issued in a run with converged ionic\n" 596 "positions; ionic (and lattice) minimization are bypassed by the vibrations module.\n" 597 "\n" 598 "Any number of the following subcommands and their arguments may follow:\n" 599 "+ dr <dr>: perturbation amplitude in bohrs for force matrix calculation (default: 0.01).\n" 600 "+ centralDiff yes|no: use a central difference formula for the second derivative\n" 601 " to achieve higher accuracy at twice the cost (default: no)\n" 602 "+ useConstraints yes|no: restrict modes of motion as specified by move flags\n" 603 " and constraints in the ion command (default: no)\n" 604 "+ translationSym yes|no: whether to assume overall translation symmetry (default yes).\n" 605 " Can be turned off to get vibrational levels in an external potential.\n" 606 "+ rotationSym yes|no: project out rotational modes (default no). Improves reliability for\n" 607 " molecular calculations. Valid only for geometries with an unambiguous center of mass.\n" 608 "+ omegaMin <omegaMin>: frequency cutoff (in Eh) for free energy calculation (default: 2e-4)\n" 609 "+ T <T>: temperature (in Kelvin) for free energy calculation (default: 298)\n" 610 "+ omegaResolution <omegaResolution>: resolution for detecting and reporting degeneracies\n" 611 " in modes (default: 1e-4). Does not affect free energies and all modes are still printed.\n" 612 "\n" 613 "Note that for a periodic system with k-points, wave functions may be incompatible\n" 614 "with and without the vibrations command due to symmetry-breaking by the perturbations.\n" 615 "To avoid this, perform electronic optimization (without initial-state) in the\n" 616 "vibration calculation itself, or consider using the phonon code instead." 617 ; 618 forbid("fix-electron-density"); 619 forbid("fix-electron-potential"); 620 } 621 processCommandVibrations622 void process(ParamList& pl, Everything& e) 623 { e.vibrations = std::make_shared<Vibrations>(); 624 while(true) 625 { VibrationsMember key; 626 pl.get(key, VM_Delim, vibMap, "key"); 627 switch(key) 628 { case VM_dr: pl.get(e.vibrations->dr, 0.01, "dr", true); break; 629 case VM_centralDiff: pl.get(e.vibrations->centralDiff, false, boolMap, "centralDiff", true); break; 630 case VM_useConstraints: pl.get(e.vibrations->useConstraints, false, boolMap, "useConstraints", true); break; 631 case VM_translationSym: pl.get(e.vibrations->translationSym, true, boolMap, "translationSym", true); break; 632 case VM_rotationSym: pl.get(e.vibrations->rotationSym, false, boolMap, "rotationSym", true); break; 633 case VM_omegaMin: pl.get(e.vibrations->omegaMin, 2e-4, "omegaMin", true); break; 634 case VM_T: pl.get(e.vibrations->T, 298., "T", true); e.vibrations->T *= Kelvin; break; 635 case VM_omegaResolution: pl.get(e.vibrations->omegaResolution, 1e-4, "omegaResolution", true); break; 636 case VM_Delim: return; //end of input 637 } 638 } 639 640 } 641 printStatusCommandVibrations642 void printStatus(Everything& e, int iRep) 643 { logPrintf("\\\n\tdr %g", e.vibrations->dr); 644 logPrintf("\\\n\tcentralDiff %s", boolMap.getString(e.vibrations->centralDiff)); 645 logPrintf("\\\n\tuseConstraints %s", boolMap.getString(e.vibrations->useConstraints)); 646 logPrintf("\\\n\ttranslationSym %s", boolMap.getString(e.vibrations->translationSym)); 647 logPrintf("\\\n\trotationSym %s", boolMap.getString(e.vibrations->rotationSym)); 648 logPrintf("\\\n\tomegaMin %g", e.vibrations->omegaMin); 649 logPrintf("\\\n\tT %g", e.vibrations->T/Kelvin); 650 logPrintf("\\\n\tomegaResolution %g", e.vibrations->omegaResolution); 651 } 652 } 653 commandVibrations; 654 655 //------------------------------------------------------------------------------------------------- 656 657 struct CommandSlabEpsilon : public Command 658 { CommandSlabEpsilonCommandSlabEpsilon659 CommandSlabEpsilon() : Command("slab-epsilon", "jdftx/Output") 660 { 661 format = "<DtotFile> <sigma> [<Ex>=0] [<Ey>=0] [<Ez>=0]"; 662 comments = 663 "Calculate dielectric function of a slab given the electrostatic potential\n" 664 "output from another calculation on same system with a different electric field.\n" 665 "+ <DtotFile> contains the electrostatic potential from the other calculation\n" 666 "+ <sigma> in bohrs specifies gaussian smoothing in output\n" 667 "+ optional <Ex>,<Ey>,Ez> specify the electric-field applied\n" 668 " in the calculation that generated <DtotFile>."; 669 670 require("coulomb-truncation-embed"); 671 } 672 processCommandSlabEpsilon673 void process(ParamList& pl, Everything& e) 674 { if(e.coulombParams.geometry != CoulombParams::Slab) 675 throw string("coulomb-interaction must be in Slab mode"); 676 e.dump.slabEpsilon = std::make_shared<SlabEpsilon>(); 677 SlabEpsilon& se = *(e.dump.slabEpsilon); 678 pl.get(se.dtotFname, string(), "DtotFile", true); 679 pl.get(se.sigma, 0., "sigma", true); 680 pl.get(se.Efield[0], 0., "Ex"); 681 pl.get(se.Efield[1], 0., "Ey"); 682 pl.get(se.Efield[2], 0., "Ez"); 683 if(!(se.Efield-e.coulombParams.Efield).length_squared()) 684 throw(string("Applied electric fields in reference and present calculations are equal")); 685 e.dump.insert(std::make_pair(DumpFreq_End, DumpSlabEpsilon)); //dump at end by default 686 } 687 printStatusCommandSlabEpsilon688 void printStatus(Everything& e, int iRep) 689 { const SlabEpsilon& se = *(e.dump.slabEpsilon); 690 logPrintf("%s %lg %lg %lg %lg", se.dtotFname.c_str(), se.sigma, se.Efield[0], se.Efield[1], se.Efield[2]); 691 } 692 } 693 commandSlabEpsilon; 694 695 //------------------------------------------------------------------------------------------------- 696 697 struct CommandBulkEpsilon : public Command 698 { CommandBulkEpsilonCommandBulkEpsilon699 CommandBulkEpsilon() : Command("bulk-epsilon", "jdftx/Output") 700 { 701 format = "<DtotFile> [<Ex>=0] [<Ey>=0] [<Ez>=0]"; 702 comments = 703 "Calculate dielectric constant of a bulk material given the electrostatic potential\n" 704 "output from another calculation on same system with a different electric field.\n" 705 "+ <DtotFile> contains the electrostatic potential from the other calculation\n" 706 "+ optional <Ex>,<Ey>,Ez> specify the electric-field applied\n" 707 " in the calculation that generated <DtotFile>.\n" 708 "It is recommended to apply field only to one reciprocal lattice direction,\n" 709 "and use a supercell of the material along that direction."; 710 711 require("electric-field"); 712 } 713 processCommandBulkEpsilon714 void process(ParamList& pl, Everything& e) 715 { if(e.coulombParams.geometry != CoulombParams::Periodic) 716 throw string("coulomb-interaction must be in Periodic mode"); 717 e.dump.bulkEpsilon = std::make_shared<BulkEpsilon>(); 718 BulkEpsilon& be = *(e.dump.bulkEpsilon); 719 pl.get(be.dtotFname, string(), "DtotFile", true); 720 pl.get(be.Efield[0], 0., "Ex"); 721 pl.get(be.Efield[1], 0., "Ey"); 722 pl.get(be.Efield[2], 0., "Ez"); 723 if(!(be.Efield-e.coulombParams.Efield).length_squared()) 724 throw(string("Applied electric fields in reference and present calculations are equal")); 725 e.dump.insert(std::make_pair(DumpFreq_End, DumpBulkEpsilon)); //dump at end by default 726 } 727 printStatusCommandBulkEpsilon728 void printStatus(Everything& e, int iRep) 729 { const BulkEpsilon& be = *(e.dump.bulkEpsilon); 730 logPrintf("%s %lg %lg %lg", be.dtotFname.c_str(), be.Efield[0], be.Efield[1], be.Efield[2]); 731 } 732 } 733 commandBulkEpsilon; 734 735 //------------------------------------------------------------------------------------------------- 736 737 extern EnumStringMap<int> truncationDirMap; 738 739 struct CommandChargedDefectCorrection : public Command 740 { CommandChargedDefectCorrectionCommandChargedDefectCorrection741 CommandChargedDefectCorrection() : Command("charged-defect-correction", "jdftx/Output") 742 { 743 format = "[Slab <dir>=100|010|001] <DtotFile> <bulkEps>|<slabEpsFile> <rMin> <rSigma>"; 744 comments = 745 "Calculate energy correction for bulk or surface charged defects \\cite ElectrostaticPotential\n" 746 "The correction is calculated assuming the defects to be model\n" 747 "charges specified using command charged-defect.\n" 748 "\n" 749 "By default, the defect is assumed bulk for coulomb-interaction Periodic\n" 750 "and surface for coulomb-interaction Slab. However, for the Periodic case,\n" 751 "the optional [Slab <dir>] overrides this to calculate surface defects\n" 752 "without truncated Coulomb potentials. Note that coulomb-truncation-embed\n" 753 "must be specified when using truncated coulomb potentials in Slab mode.\n" 754 "In Periodic mode, the correction assumes a slab centered at the origin\n" 755 "(i.e. analogous to using xCenter 0 0 0 in the truncated mode).\n" 756 "\n" 757 "<DtotFile> contains the electrostatic potential from a reference\n" 758 "neutral calculation with similar geometry (lattice vectors and grid\n" 759 "must match exactly).\n" 760 "\n" 761 "For bulk defect calculations, <bulkEps> is the bulk dielectric constant.\n" 762 "\n" 763 "For surface defect calculations, <slabEpsFile> specifies a dielectric\n" 764 "profile calculated using command slab-epsilon in a similar geometry\n" 765 "(the number of points along the slab normal direction must match exactly).\n" 766 "Optionally, the <slabEpsFile> may contain an additional column for the\n" 767 "in-plane response (which is not computed by command slab-epsilon),\n" 768 "in which case an anisotropic dielectric model is used.\n" 769 "\n" 770 "<rMin> specifies the distance away from the defect center to use in\n" 771 "the determination of the alignment potential, with rSigma specifying an\n" 772 "error function turn-on distance. The code wil generate a text file with\n" 773 "the spherically averaged model and DFT electrostatic potentials, which\n" 774 "can be used to check the calculated alignment and refine rMin and rSigma."; 775 776 require("latt-scale"); 777 require("coords-type"); 778 require("coulomb-interaction"); 779 } 780 processCommandChargedDefectCorrection781 void process(ParamList& pl, Everything& e) 782 { e.dump.chargedDefect = std::make_shared<ChargedDefect>(); 783 ChargedDefect& cd = *(e.dump.chargedDefect); 784 //Check for geometry override: 785 cd.geometry = e.coulombParams.geometry; 786 cd.iDir = e.coulombParams.iDir; 787 string slabSpec; 788 pl.get(slabSpec, string(), "Slab|DtotFile", true); 789 if(slabSpec == "Slab") 790 { if(cd.geometry != CoulombParams::Periodic) 791 throw string("Slab geometry override should only be specified for coulomb-interaction Periodic"); 792 cd.geometry = CoulombParams::Slab; 793 pl.get(cd.iDir, 0, truncationDirMap, "dir", true); 794 } 795 else pl.rewind(); 796 //Ref potential: 797 pl.get(cd.dtotFname, string(), "DtotFile", true); 798 //Dielectric function (and add citation depending on geometry): 799 switch(cd.geometry) 800 { case CoulombParams::Periodic: 801 { pl.get(cd.bulkEps, 1., "bulkEps", true); 802 Citations::add("Correction scheme for charged bulk defects", 803 "C Freysoldt, J Neugebauer and C. van de Walle, Phys. Rev. Lett. 102, 016402 (2009)"); 804 break; 805 } 806 case CoulombParams::Slab: 807 { pl.get(cd.slabEpsFname, string(), "slabEpsFile", true); 808 Citations::add("Correction scheme for charged surface defects", 809 "H Komsa and A Pasquarello, Alfredo, Phys. Rev. Lett. 110, 095505 (2013)"); 810 break; 811 } 812 default: throw string("coulomb-interaction must be either Slab or Periodic"); 813 } 814 //Alignment potential ranges: 815 pl.get(cd.rMin, 0., "rMin", true); 816 pl.get(cd.rSigma, 0., "rSigma", true); 817 e.dump.insert(std::make_pair(DumpFreq_End, DumpChargedDefect)); //dump at end by default 818 } 819 printStatusCommandChargedDefectCorrection820 void printStatus(Everything& e, int iRep) 821 { const ChargedDefect& cd = *(e.dump.chargedDefect); 822 if(cd.geometry != e.coulombParams.geometry) 823 logPrintf("Slab %s ", truncationDirMap.getString(cd.iDir)); 824 logPrintf("%s ", cd.dtotFname.c_str()); 825 switch(cd.geometry) 826 { case CoulombParams::Periodic: logPrintf("%lg", cd.bulkEps); break; 827 case CoulombParams::Slab: logPrintf("%s", cd.slabEpsFname.c_str()); break; 828 default:; //should never be encountered 829 } 830 logPrintf(" %lg %lg", cd.rMin, cd.rSigma); 831 } 832 } 833 CommandChargedDefectCorrection; 834 835 struct CommandChargedDefect : public Command 836 { CommandChargedDefectCommandChargedDefect837 CommandChargedDefect() : Command("charged-defect", "jdftx/Output") 838 { 839 format = "<x0> <x1> <x2> <q> <sigma>"; 840 comments = 841 "Specify model charge distribution(s) for charged-defect-correction.\n" 842 "Issue the command once for each defect in the calculation cell.\n" 843 "The defect charge density will be modeled by a Gaussian of width\n" 844 "<sigma> with net electron count <q> located at <x0>,<x1>,<x2>\n" 845 "(in the coordinate system specified by coords-type). The Gaussian must\n" 846 "be resolvable on the plane-wave grid; recommend rSigma >= 1 bohr."; 847 848 require("charged-defect-correction"); 849 allowMultiple = true; 850 } 851 processCommandChargedDefect852 void process(ParamList& pl, Everything& e) 853 { ChargedDefect::Center cdc; 854 //Position: 855 vector3<> pos; 856 pl.get(pos[0], 0., "x0", true); 857 pl.get(pos[1], 0., "x1", true); 858 pl.get(pos[2], 0., "x2", true); 859 cdc.pos = (e.iInfo.coordsType==CoordsLattice) ? pos : inv(e.gInfo.R)*pos; 860 //Charge and width: 861 pl.get(cdc.q, 0., "q", true); 862 pl.get(cdc.sigma, 0., "sigma", true); 863 e.dump.chargedDefect->center.push_back(cdc); 864 } 865 printStatusCommandChargedDefect866 void printStatus(Everything& e, int iRep) 867 { const ChargedDefect::Center& cdc = e.dump.chargedDefect->center[iRep]; 868 vector3<> pos = (e.iInfo.coordsType==CoordsLattice) ? cdc.pos : e.gInfo.R*cdc.pos; 869 logPrintf("%lg %lg %lg %+lg %lg", pos[0], pos[1], pos[2], cdc.q, cdc.sigma); 870 } 871 } 872 commandChargedDefect; 873 874 struct CommandPotentialSubtraction : public Command 875 { CommandPotentialSubtractionCommandPotentialSubtraction876 CommandPotentialSubtraction() : Command("potential-subtraction", "jdftx/Output") 877 { format = "<subtract>=yes|no"; 878 comments = 879 "Whether to subtract neutral atom potentials in dumped potentials (Dtot and Dvac).\n" 880 "This subtraction produces much smoother potentials and is enabled by default \\cite ElectrostaticPotential."; 881 } 882 processCommandPotentialSubtraction883 void process(ParamList& pl, Everything& e) 884 { pl.get(e.dump.potentialSubtraction, true, boolMap, "subtract"); 885 } 886 printStatusCommandPotentialSubtraction887 void printStatus(Everything& e, int iRep) 888 { logPrintf("%s", boolMap.getString(e.dump.potentialSubtraction)); 889 } 890 } 891 commandPotentialSubtraction; 892 893 894 struct CommandBandUnfold : public Command 895 { CommandBandUnfoldCommandBandUnfold896 CommandBandUnfold() : Command("band-unfold", "jdftx/Output") 897 { 898 format = " \\\n\t<M00> <M01> <M02> \\\n\t<M10> <M11> <M12> \\\n\t<M20> <M21> <M22>"; 899 comments = 900 "Unfold band structure from a supercell calculation to a unit cell\n" 901 "with lattice vectors Runit, defined by the integer matrix M such\n" 902 "that current lattice vectors R = Runit * M."; 903 } 904 processCommandBandUnfold905 void process(ParamList& pl, Everything& e) 906 { matrix3<int>& M = e.dump.Munfold; 907 for(int j=0; j<3; j++) for(int k=0; k<3; k++) 908 { ostringstream oss; oss << "s" << j << k; 909 pl.get(M(j,k), 0, oss.str(), true); 910 } 911 e.dump.insert(std::make_pair(DumpFreq_End, DumpBandUnfold)); 912 } 913 printStatusCommandBandUnfold914 void printStatus(Everything& e, int iRep) 915 { for (int j=0; j < 3; j++) 916 { logPrintf(" \\\n\t"); 917 for(int k=0; k<3; k++) logPrintf("%d ",e.dump.Munfold(j,k)); 918 } 919 } 920 } 921 commandBandUnfold; 922 923 924 enum BGWparamsMember 925 { BGWpm_nBandsDense, 926 BGWpm_blockSize, 927 BGWpm_clusterSize, 928 BGWpm_EcutChiFluid, 929 BGWpm_elecOnly, 930 BGWpm_q0, 931 BGWpm_freqReMax_eV, 932 BGWpm_freqReStep_eV, 933 BGWpm_freqBroaden_eV, 934 BGWpm_freqNimag, 935 BGWpm_freqPlasma, 936 BGWpm_Delim 937 }; 938 EnumStringMap<BGWparamsMember> bgwpmMap 939 ( BGWpm_nBandsDense, "nBandsDense", 940 BGWpm_blockSize, "blockSize", 941 BGWpm_clusterSize, "clusterSize", 942 BGWpm_EcutChiFluid, "EcutChiFluid", 943 BGWpm_elecOnly, "elecOnly", 944 BGWpm_q0, "q0", 945 BGWpm_freqReMax_eV, "freqReMax_eV", 946 BGWpm_freqReStep_eV, "freqReStep_eV", 947 BGWpm_freqBroaden_eV, "freqBroaden_eV", 948 BGWpm_freqNimag, "freqNimag", 949 BGWpm_freqPlasma, "freqPlasma" 950 ); 951 EnumStringMap<BGWparamsMember> bgwpmDescMap 952 ( BGWpm_nBandsDense, "If non-zero, use a dense ScaLAPACK solver to calculate more bands", 953 BGWpm_blockSize, "Block size for ScaLAPACK diagonalization (default: 32)", 954 BGWpm_clusterSize, "Maximum eigenvalue cluster size to allocate extra ScaLAPACK workspace for (default: 10)", 955 BGWpm_EcutChiFluid, "KE cutoff in hartrees for fluid polarizability output (default: 0; set non-zero to enable)", 956 BGWpm_elecOnly, "Whether fluid polarizability output should only include electronic response (default: true)", 957 BGWpm_q0, "Zero wavevector replacement to be used for polarizability output (default: (0,0,0))", 958 BGWpm_freqReMax_eV, "Maximum real frequency in eV (default: 30.)", 959 BGWpm_freqReStep_eV, "Real frequency grid spacing in eV (default: 1.)", 960 BGWpm_freqBroaden_eV, "Broadening (imaginary part) of real frequency grid in eV (default: 0.1)", 961 BGWpm_freqNimag, "Number of imaginary frequencies (default: 25)", 962 BGWpm_freqPlasma, "Plasma frequency in Hartrees used in GW imaginary frequency grid (default: 1.), set to zero for RPA frequency grid" 963 ); 964 965 struct CommandBGWparams : public Command 966 { CommandBGWparamsCommandBGWparams967 CommandBGWparams() : Command("bgw-params", "jdftx/Output") 968 { 969 format = "<key1> <value1> <key2> <value2> ..."; 970 comments = "Control BGW output. Possible keys and value types are:" 971 + addDescriptions(bgwpmMap.optionList(), linkDescription(bgwpmMap, bgwpmDescMap)) 972 + "\n\nAny number of these key-value pairs may be specified in any order."; 973 } 974 processCommandBGWparams975 void process(ParamList& pl, Everything& e) 976 { e.dump.bgwParams = std::make_shared<BGWparams>(); 977 BGWparams& bgwp = *(e.dump.bgwParams); 978 while(true) 979 { BGWparamsMember key; 980 pl.get(key, BGWpm_Delim, bgwpmMap, "key"); 981 #define READ_AND_CHECK(param, op, val) \ 982 case BGWpm_##param: \ 983 pl.get(bgwp.param, val, #param, true); \ 984 if(!(bgwp.param op val)) throw string(#param " must be " #op " " #val); \ 985 break; 986 switch(key) 987 { READ_AND_CHECK(nBandsDense, >=, 0) 988 READ_AND_CHECK(blockSize, >, 0) 989 READ_AND_CHECK(clusterSize, >, 0) 990 READ_AND_CHECK(EcutChiFluid, >=, 0.) 991 case BGWpm_elecOnly: 992 pl.get(bgwp.elecOnly, true, boolMap, "elecOnly", true); 993 break; 994 case BGWpm_q0: 995 for(int dir=0; dir<3; dir++) 996 pl.get(bgwp.q0[dir], 0., "q0", true); 997 break; 998 READ_AND_CHECK(freqReMax_eV, >, 0.) 999 READ_AND_CHECK(freqReStep_eV, >, 0.) 1000 READ_AND_CHECK(freqBroaden_eV, >, 0.) 1001 READ_AND_CHECK(freqNimag, >, 0) 1002 READ_AND_CHECK(freqPlasma, >=, 0.) 1003 case BGWpm_Delim: return; //end of input 1004 } 1005 #undef READ_AND_CHECK 1006 } 1007 } 1008 printStatusCommandBGWparams1009 void printStatus(Everything& e, int iRep) 1010 { assert(e.dump.bgwParams); 1011 const BGWparams& bgwp = *(e.dump.bgwParams); 1012 #define PRINT(param, format) logPrintf(" \\\n\t" #param " " format, bgwp.param); 1013 PRINT(nBandsDense, "%d") 1014 PRINT(blockSize, "%d") 1015 PRINT(clusterSize, "%d") 1016 PRINT(EcutChiFluid, "%lg") 1017 logPrintf(" \\\n\telecOnly %s", boolMap.getString(bgwp.elecOnly)); 1018 logPrintf(" \\\n\tq0 %lg %lg %lg", bgwp.q0[0], bgwp.q0[1], bgwp.q0[2]); 1019 PRINT(freqReMax_eV, "%lg") 1020 PRINT(freqReStep_eV, "%lg") 1021 PRINT(freqBroaden_eV, "%lg") 1022 PRINT(freqNimag, "%d") 1023 PRINT(freqPlasma, "%lg") 1024 #undef PRINT 1025 } 1026 } 1027 commandBGWparams; 1028