1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth Weaver
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 <core/Units.h>
23
24 EnumStringMap<FluidType> fluidTypeMap
25 ( FluidNone, "None",
26 FluidLinearPCM, "LinearPCM",
27 FluidNonlinearPCM, "NonlinearPCM",
28 FluidSaLSA, "SaLSA",
29 FluidClassicalDFT, "ClassicalDFT"
30 );
31
32 struct CommandFluid : public Command
33 {
CommandFluidCommandFluid34 CommandFluid() : Command("fluid", "jdftx/Fluid/Parameters")
35 {
36 format = "[<type>=None] [<Temperature>=298K] [<Pressure>=1.01325bar]";
37 comments = "Perform joint density functional theory with fluid of <type>:\n"
38 "\n+ None:\n\n"
39 " Standard vacuum DFT calculation with no solvation model.\n"
40 "\n+ LinearPCM: \\cite NonlinearPCM \\cite CANDLE \\cite PCM-SCCS\n\n"
41 " Use a solvation model that includes linear dielectric (and/or ionic)\n"
42 " response. Select a specific linear solvation model using pcm-variant.\n"
43 "\n+ NonlinearPCM: \\cite NonlinearPCM \\cite CavityWDA\n\n"
44 " Use a solvation model that includes nonlinear dielectric (and/or ionic)\n"
45 " response, and accounts for dielectric saturation effects.\n"
46 " Select a specific nonlinear solvation model using pcm-variant.\n"
47 "\n+ SaLSA: \\cite SaLSA\n\n"
48 " Use the non-empirical nonlocal-response solvation model based on the\n"
49 " Spherically-averaged Liquid Susceptibility Ansatz.\n"
50 "\n+ ClassicalDFT: \\cite PolarizableCDFT \\cite RigidCDFT \\cite BondedVoids\n\n"
51 " Full joint density-functional theory with a classical density-functional\n"
52 " description of the solvent. See fluid-solvent, fluid-cation, fluid-anion\n"
53 " and related commands for controlling the classical density-functional theory.\n"
54 "\n"
55 "Optionally adjust the fluid <Temperature> (in Kelvin) and <Pressure> (in bars).";
56 hasDefault = true;
57 require("coulomb-interaction");
58 }
59
processCommandFluid60 void process(ParamList& pl, Everything& e)
61 { FluidSolverParams& fsp = e.eVars.fluidParams;
62 pl.get(fsp.fluidType, FluidNone, fluidTypeMap, "type");
63 if((e.coulombParams.geometry != CoulombParams::Periodic) && (fsp.fluidType != FluidNone))
64 e.coulombParams.embedFluidMode = true; //require embedding in fluid mode (periodic Coulomb kernels in larger box)
65 pl.get(fsp.T, 298., "Temperature"); fsp.T *= Kelvin; //convert to atomic units
66 pl.get(fsp.P, 1.01325, "Pressure"); fsp.P *= Bar; //convert to atomic units
67 }
68
printStatusCommandFluid69 void printStatus(Everything& e, int iRep)
70 { const FluidSolverParams& fsp = e.eVars.fluidParams;
71 logPrintf("%s", fluidTypeMap.getString(fsp.fluidType));
72 if(fsp.fluidType != FluidNone)
73 logPrintf(" %lf %lf", fsp.T/Kelvin, fsp.P/Bar);
74 }
75 }
76 commandFluid;
77
78
79 struct CommandFluidGummelLoop : public Command
80 {
CommandFluidGummelLoopCommandFluidGummelLoop81 CommandFluidGummelLoop() : Command("fluid-gummel-loop", "jdftx/Fluid/Optimization")
82 {
83 format = "[<maxIterations>=10] [<Atol>=1e-5]";
84 comments =
85 "Settings for the fluid <--> electron self-consistency loop:\n"
86 "+ <maxIterations>: Max number of electron and fluid minimization pairs\n"
87 "+ <Atol>: Free energy convergence criterion for this outer loop.\n"
88 "Use fluid-solve-frequency to control whether such a loop is used at all.";
89 hasDefault = true;
90 }
91
processCommandFluidGummelLoop92 void process(ParamList& pl, Everything& e)
93 { pl.get(e.cntrl.fluidGummel_nIterations, 10, "maxIterations");
94 pl.get(e.cntrl.fluidGummel_Atol, 1e-5, "Atol");
95 }
96
printStatusCommandFluidGummelLoop97 void printStatus(Everything& e, int iRep)
98 { logPrintf("%d %le", e.cntrl.fluidGummel_nIterations, e.cntrl.fluidGummel_Atol);
99 }
100 }
101 commandFluidGummelLoop;
102
103
104 EnumStringMap<FluidSolveFrequency> fluidSolveFreqMap
105 ( FluidFreqInner, "Inner",
106 FluidFreqGummel, "Gummel",
107 FluidFreqDefault, "Default"
108 );
109 EnumStringMap<FluidSolveFrequency> fluidSolveFreqDescMap
110 ( FluidFreqInner, "Solve fluid every electronic step",
111 FluidFreqGummel, "Alternately minimize fluid and electrons (fluid-gummel-loop)",
112 FluidFreqDefault, "Decide based on fluid type (Inner for LinearPCM/SaLSA, Gummel for rest)"
113 );
114
115 struct CommandFluidSolveFrequency : public Command
116 {
CommandFluidSolveFrequencyCommandFluidSolveFrequency117 CommandFluidSolveFrequency() : Command("fluid-solve-frequency", "jdftx/Fluid/Optimization")
118 {
119 format = "<freq>=" + fluidSolveFreqMap.optionList();
120 comments = "Select how often to optimize fluid state:"
121 + addDescriptions(fluidSolveFreqMap.optionList(), linkDescription(fluidSolveFreqMap, fluidSolveFreqDescMap));
122
123 require("fluid");
124 }
125
processCommandFluidSolveFrequency126 void process(ParamList& pl, Everything& e)
127 { FluidSolverParams& fsp = e.eVars.fluidParams;
128 pl.get(fsp.solveFrequency, FluidFreqDefault, fluidSolveFreqMap, "freq", true);
129 //Check for cases that don't support Gummel loop:
130 if(fsp.solveFrequency==FluidFreqGummel)
131 { if(fsp.fluidType==FluidLinearPCM) throw string("Fluid type 'LinearPCM' does not support fluid-solve-frequency Gummel");
132 if(fsp.fluidType==FluidSaLSA) throw string("Fluid type 'SaLSA' does not support fluid-solve-frequency Gummel");
133 }
134 }
135
printStatusCommandFluidSolveFrequency136 void printStatus(Everything& e, int iRep)
137 { const FluidSolverParams& fsp = e.eVars.fluidParams;
138 logPrintf("%s", fluidSolveFreqMap.getString(fsp.solveFrequency));
139 }
140 }
141 commandFluidSolveFrequency;
142
143
144 struct CommandFluidInitialState : public Command
145 {
CommandFluidInitialStateCommandFluidInitialState146 CommandFluidInitialState() : Command("fluid-initial-state", "jdftx/Initialization")
147 {
148 format = "<filename>";
149 comments = "Read initial state of a fluid (compatible with *.fluidState from dump End State)";
150
151 forbid("initial-state");
152 }
153
processCommandFluidInitialState154 void process(ParamList& pl, Everything& e)
155 { pl.get(e.eVars.fluidInitialStateFilename, string(), "filename", true);
156 }
157
printStatusCommandFluidInitialState158 void printStatus(Everything& e, int iRep)
159 { logPrintf("%s", e.eVars.fluidInitialStateFilename.c_str());
160 }
161 }
162 commandFluidInitialState;
163
164
165 struct CommandFluidVdwScale : public Command
166 {
CommandFluidVdwScaleCommandFluidVdwScale167 CommandFluidVdwScale() : Command("fluid-vdwScale", "jdftx/Fluid/Parameters")
168 {
169 format = "<scale=0.75>";
170 comments = "Scale van der Waals interactions between fluid and explicit system by a constant factor <scale>.\n\n"
171 "Default is fluid specific and ranges between 0.4 to 1.3.\n"
172 "Set to 0 to use the prefactor corresponding to fluid exchange-correlation.";
173 require("fluid-solvent");
174 }
175
processCommandFluidVdwScale176 void process(ParamList& pl, Everything& e)
177 {
178 pl.get(e.eVars.fluidParams.vdwScale, 0.75, "scale");
179 }
180
printStatusCommandFluidVdwScale181 void printStatus(Everything& e, int iRep)
182 {
183 logPrintf(" %lg", e.eVars.fluidParams.vdwScale);
184 }
185 }
186 commandFluidVDWscale;
187
188
189 EnumStringMap<FluidComponent::Name> solventMap
190 ( FluidComponent::H2O, "H2O",
191 FluidComponent::CHCl3, "CHCl3",
192 FluidComponent::CCl4, "CCl4",
193 FluidComponent::CH3CN, "CH3CN",
194 FluidComponent::DMC, "DMC",
195 FluidComponent::EC, "EC",
196 FluidComponent::PC, "PC",
197 FluidComponent::DMF, "DMF",
198 FluidComponent::THF, "THF",
199 FluidComponent::DMSO, "DMSO",
200 FluidComponent::CH2Cl2, "CH2Cl2",
201 FluidComponent::Ethanol, "Ethanol",
202 FluidComponent::Methanol, "Methanol",
203 FluidComponent::Octanol, "Octanol",
204 FluidComponent::EthylEther, "EthylEther",
205 FluidComponent::Chlorobenzene, "Chlorobenzene",
206 FluidComponent::Isobutanol, "Isobutanol",
207 FluidComponent::CarbonDisulfide, "CarbonDisulfide",
208 FluidComponent::Glyme, "Glyme",
209 FluidComponent::EthyleneGlycol, "EthyleneGlycol"
210 );
211 EnumStringMap<FluidComponent::Name> cationMap
212 ( FluidComponent::Sodium, "Na+",
213 // FluidComponent::HydratedSodium, "Na(6H2O)+",
214 FluidComponent::Potassium, "K+"
215 // FluidComponent::HydratedPotassium, "K(6H2O)+",
216 // FluidComponent::Hydronium, "H3O+",
217 // FluidComponent::HydratedHydronium, "H3O(4H2O)+"
218 );
219 EnumStringMap<FluidComponent::Name> anionMap
220 ( FluidComponent::Chloride, "Cl-",
221 FluidComponent::Fluoride, "F-",
222 FluidComponent::Perchlorate, "ClO4-"
223 // FluidComponent::Hydroxide, "OH-",
224 // FluidComponent::HydratedHydroxide, "OH(4H2O)-"
225 );
226 EnumStringMap<FluidComponent::Functional> functionalMap
227 ( FluidComponent::ScalarEOS, "ScalarEOS",
228 FluidComponent::FittedCorrelations, "FittedCorrelations",
229 FluidComponent::BondedVoids, "BondedVoids",
230 FluidComponent::MeanFieldLJ, "MeanFieldLJ"
231 );
232
233 EnumStringMap<FluidComponent::TranslationMode> translationModeMap
234 ( FluidComponent::ConstantSpline, "ConstantSpline",
235 FluidComponent::LinearSpline, "LinearSpline",
236 FluidComponent::Fourier, "Fourier"
237 );
238
239 EnumStringMap<FluidComponent::Representation> representationMap
240 ( FluidComponent::MuEps, "MuEps",
241 FluidComponent::Pomega, "Pomega",
242 FluidComponent::PsiAlpha, "PsiAlpha"
243 );
244
245 const EnumStringMap<S2quadType>& s2quadTypeMap = S2quadTypeMap;
246
247 enum FluidComponentMember
248 { FCM_epsBulk, //!< bulk dielectric constant
249 FCM_pMol, //!< dipole moment of each molecule in e-bohr
250 FCM_epsInf, //!< optical-frequency dielectric constant
251 FCM_Pvap, //!< vapor pressure in Eh/bohr^3
252 FCM_sigmaBulk, //!< bulk surface tension in Eh/bohr^2
253 FCM_Rvdw, //!< effective van der Waals radius of the fluid (derived from equation of state) in bohrs
254 FCM_Res, //!< electrostatic radius of solvent (derived from nonlocal response) in bohrs
255 //Extras for frequency dependence:
256 FCM_tauNuc, //!< nuclear motion damping time (in Eh^-1 atomic units): rotational for solvents, translational for ions
257 FCM_poleEl, //!< electronic response poles (Drude-Lorentz model)
258 //Extras for ClassicalDFT:
259 FCM_epsLJ, //!< Lennard-Jones well depth for Mean-Field LJ excess functional
260 FCM_representation, //!< ideal gas representation
261 FCM_s2quadType, //!< orientation quadrature type
262 FCM_quad_nBeta, //!< number of beta samples for Euler quadrature
263 FCM_quad_nAlpha, //!< number of alpha samples for Euler quadrature
264 FCM_quad_nGamma, //!< number of gamma samples for Euler quadrature
265 FCM_translationMode, //!< translation operator type
266 FCM_Nnorm, //!< unit cell molecule count constraint
267 //Delimiter used in parsing:
268 FCM_Delim
269 };
270 EnumStringMap<FluidComponentMember> fcmMap
271 ( FCM_epsBulk, "epsBulk",
272 FCM_pMol, "pMol",
273 FCM_epsInf, "epsInf",
274 FCM_Pvap, "Pvap",
275 FCM_sigmaBulk, "sigmaBulk",
276 FCM_Rvdw, "Rvdw",
277 FCM_Res, "Res",
278 FCM_tauNuc, "tauNuc",
279 FCM_poleEl, "poleEl",
280 FCM_epsLJ, "epsLJ",
281 FCM_representation, "representation",
282 FCM_s2quadType, "s2quadType",
283 FCM_quad_nBeta, "quad_nBeta",
284 FCM_quad_nAlpha, "quad_nAlpha",
285 FCM_quad_nGamma, "quad_nGamma",
286 FCM_translationMode,"translation",
287 FCM_Nnorm, "Nnorm"
288 );
289 EnumStringMap<FluidComponentMember> fcmDescMap
290 ( FCM_epsBulk, "bulk dielectric constant",
291 FCM_pMol, "dipole moment of each molecule in e-bohr",
292 FCM_epsInf, "optical-frequency dielectric constant",
293 FCM_Pvap, "vapor pressure in Eh/bohr^3",
294 FCM_sigmaBulk, "bulk surface tension in Eh/bohr^2",
295 FCM_Rvdw, "effective van der Waals radius of the fluid (derived from equation of state) in bohrs",
296 FCM_Res, "electrostatic radius of solvent (derived from nonlocal response) in bohrs",
297 FCM_tauNuc, "nuclear motion damping time (in Eh^-1 atomic units): rotational for solvents, translational for ions",
298 FCM_poleEl, "electronic response Lorentz poles with parameters ( omega0[eV] gamma0[eV] A0 ). [specify multiple times for several poles, with A0 adding up to 1]",
299 FCM_epsLJ, "Lennard-Jones well depth for Mean-Field LJ excess functional",
300 FCM_representation, "ideal gas representation: " + addDescriptions(representationMap.optionList(), nullDescription, "\n - "),
301 FCM_s2quadType, "orientation quadrature type:" + addDescriptions(s2quadTypeMap.optionList(), nullDescription, "\n - "),
302 FCM_quad_nBeta, "number of beta samples for Euler quadrature",
303 FCM_quad_nAlpha, "number of alpha samples for Euler quadrature",
304 FCM_quad_nGamma, "number of gamma samples for Euler quadrature",
305 FCM_translationMode, "translation operator type: " + addDescriptions(translationModeMap.optionList(), nullDescription, "\n - "),
306 FCM_Nnorm, "unit cell molecule count constraint"
307 );
308
309 //Abstract base class for fluid-solvent fluid-cation and fluid-anion
310 struct CommandFluidComponent : public Command
311 {
312 private:
313 const EnumStringMap<FluidComponent::Name>& nameMap;
314 FluidComponent::Name defaultName;
315 FluidComponent::Functional defaultFunctional;
316 bool defaultEnabled;
317
318 protected:
CommandFluidComponentCommandFluidComponent319 CommandFluidComponent(string suffix, const EnumStringMap<FluidComponent::Name>& nameMap, FluidComponent::Name defaultName, FluidComponent::Functional defaultFunctional, bool defaultEnabled)
320 : Command("fluid-"+suffix, "jdftx/Fluid/Constituents"),
321 nameMap(nameMap), defaultName(defaultName), defaultFunctional(defaultFunctional), defaultEnabled(defaultEnabled)
322 {
323 format = (defaultEnabled ? ("[<name>=" + string(nameMap.getString(defaultName)) +"] [<concentration>=bulk]") : "<name> <concentration>")
324 + " [<functional>=" + string(functionalMap.getString(defaultFunctional)) + "]"
325 " [<key1> <value1>] ...";
326 comments = "Add " + suffix + " component to fluid, where <name> may be one of:"
327 + addDescriptions(nameMap.optionList(), nullDescription) + "\n\n"
328 + string(defaultEnabled
329 ? "The concentration may be specified explicitly in mol/liter or set to 'bulk' to use bulk fluid value."
330 : "The concentration must be specified explicitly in mol/liter.") + "\n"
331 "For classical DFT fluids, the excess functional may be one of " + functionalMap.optionList() + ".\n"
332 "\n"
333 "Optional component properties may be overridden using an arbitrary number of <key> <value> pairs.\n"
334 + (suffix=="solvent"
335 ? ("Possible keys and value types are:"
336 + addDescriptions(fcmMap.optionList(), linkDescription(fcmMap, fcmDescMap))
337 + "\n\nAny number of these key-value pairs may be specified in any order.")
338 : string("See command fluid-solvent for a description of adjustable properties."));
339
340 hasDefault = defaultEnabled;
341 allowMultiple = true;
342 }
343
344 public:
processCommandFluidComponent345 void process(ParamList& pl, Everything& e)
346 { //Read parameters:
347 FluidComponent::Name name; pl.get(name, defaultName, nameMap, "name", !defaultEnabled);
348 string keyNbulk; pl.get(keyNbulk, string("bulk"), "concentration", !defaultEnabled);
349 FluidComponent::Functional functional; pl.get(functional, defaultFunctional, functionalMap, "functional");
350 //Make component and add to list:
351 auto c = std::make_shared<FluidComponent>(name, e.eVars.fluidParams.T, functional);
352 e.eVars.fluidParams.addComponent(c);
353 //Check and optionally override concentration:
354 if((!defaultEnabled) || (keyNbulk != "bulk"))
355 { istringstream iss(keyNbulk);
356 double Nbulk = 0.; iss >> Nbulk;
357 if(iss.fail() || !iss.eof()) throw string("Conversion of parameter <concentration> failed");
358 if(Nbulk <= 0.) throw string("<concentration> must be positive");
359 c->Nbulk = Nbulk * (mol/liter);
360 }
361 //Optional properties
362 bool poleElAdded = false; //whether electronic poles have been added already by this command
363 while(true)
364 { FluidComponentMember key;
365 pl.get(key, FCM_Delim, fcmMap, "key");
366 #define READ_AND_CHECK(param, op, val) \
367 case FCM_##param: \
368 pl.get(c->param, val, #param, true); \
369 if(!(c->param op val)) throw string(#param " must be " #op " " #val); \
370 break;
371 #define READ_ENUM(param, paramDefault) \
372 case FCM_##param: \
373 pl.get(c->param, paramDefault, param##Map, #param, true); \
374 break;
375 switch(key)
376 { READ_AND_CHECK(epsBulk, >, 1.)
377 READ_AND_CHECK(pMol, >=, 0.)
378 READ_AND_CHECK(epsInf, >=, 1.)
379 READ_AND_CHECK(Pvap, >, 0.)
380 READ_AND_CHECK(sigmaBulk, >, 0.)
381 READ_AND_CHECK(Rvdw, >, 0.)
382 READ_AND_CHECK(Res, >, 0.)
383 READ_AND_CHECK(tauNuc, >, 0.)
384 case FCM_poleEl:
385 { if(!poleElAdded) c->polesEl.clear(); //remove any default poles
386 FluidComponent::PoleLD pole;
387 pl.get(pole.omega0, 0., "poleEl::omega0", true);
388 pl.get(pole.gamma0, 0., "poleEl::gamma0", true);
389 pl.get(pole.A0, 0., "poleEl::A0", true);
390 //Check:
391 if(pole.omega0 < 0.) throw string("poleEl::omega0 must be >= 0");
392 if(pole.gamma0 <= 0.) throw string("poleEl::gamma0 must be > 0");
393 pole.omega0 *= eV; //convert from eV to Eh
394 pole.gamma0 *= eV; //convert from eV to Eh
395 c->polesEl.push_back(pole);
396 poleElAdded = true;
397 break;
398 }
399 READ_AND_CHECK(epsLJ, >, 0.)
400 READ_ENUM(representation, FluidComponent::MuEps)
401 READ_ENUM(s2quadType, QuadOctahedron)
402 READ_AND_CHECK(quad_nBeta, >, 0u)
403 READ_AND_CHECK(quad_nAlpha, >=, 0u)
404 READ_AND_CHECK(quad_nGamma, >=, 0u)
405 READ_ENUM(translationMode, FluidComponent::LinearSpline)
406 READ_AND_CHECK(Nnorm, >=, 0.)
407 case FCM_Delim: return; //end of input
408 }
409 #undef READ_AND_CHECK
410 #undef READ_ENUM
411 }
412 if(poleElAdded)
413 { double A0sum = 0.;
414 for(const FluidComponent::PoleLD& pole: c->polesEl)
415 A0sum += pole.A0;
416 if(fabs(A0sum-1) > 1e-3) throw string("poleEl::A0 should add up to 1.");
417 }
418 }
419
printCommandFluidComponent420 void print(const Everything& e, const FluidComponent& c)
421 { logPrintf("%s %lg %s", nameMap.getString(c.name), c.Nbulk/(mol/liter), functionalMap.getString(c.functional));
422 #define PRINT(param) logPrintf(" \\\n\t" #param " %lg", c.param);
423 #define PRINT_UINT(param) logPrintf(" \\\n\t" #param " %u", c.param);
424 #define PRINT_ENUM(param) logPrintf(" \\\n\t" #param " %s", param##Map.getString(c.param));
425 PRINT(epsBulk)
426 PRINT(pMol)
427 PRINT(epsInf)
428 PRINT(Pvap)
429 PRINT(sigmaBulk)
430 PRINT(Rvdw)
431 PRINT(Res)
432 PRINT(tauNuc)
433 for(const FluidComponent::PoleLD& pole: c.polesEl)
434 logPrintf(" \\\n\tpoleEl %lg %lg %lg", pole.omega0/eV, pole.gamma0/eV, pole.A0);
435 if(e.eVars.fluidParams.fluidType == FluidClassicalDFT)
436 { PRINT(epsLJ)
437 PRINT_ENUM(representation)
438 PRINT_ENUM(s2quadType)
439 PRINT_UINT(quad_nBeta)
440 PRINT_UINT(quad_nAlpha)
441 PRINT_UINT(quad_nGamma)
442 PRINT_ENUM(translationMode)
443 PRINT(Nnorm)
444 }
445 #undef PRINT
446 #undef PRINT_INT
447 #undef PRINT_ENUM
448 }
449 };
450
451 struct CommandFluidSolvent : public CommandFluidComponent
452 {
CommandFluidSolventCommandFluidSolvent453 CommandFluidSolvent() : CommandFluidComponent("solvent", solventMap, FluidComponent::H2O, FluidComponent::ScalarEOS, true)
454 { require("pcm-variant"); //which in turn requires fluid
455 }
456
processCommandFluidSolvent457 void process(ParamList& pl, Everything& e)
458 { CommandFluidComponent::process(pl, e);
459 //Set PCM parameters if necessary:
460 switch(e.eVars.fluidParams.fluidType)
461 { case FluidLinearPCM:
462 case FluidNonlinearPCM:
463 case FluidSaLSA:
464 if(e.eVars.fluidParams.solvents.size()>1)
465 throw string("PCMs require exactly one solvent component - more than one specified.");
466 e.eVars.fluidParams.setPCMparams();
467 break;
468 case FluidClassicalDFT:
469 e.eVars.fluidParams.setCDFTparams();
470 break;
471 default:;
472 }
473 }
474
printStatusCommandFluidSolvent475 void printStatus(Everything& e, int iRep)
476 { print(e, *(e.eVars.fluidParams.solvents[iRep]));
477 }
478 }
479 commandFluidSolvent;
480
481 struct CommandFluidCation : public CommandFluidComponent
482 {
CommandFluidCationCommandFluidCation483 CommandFluidCation() : CommandFluidComponent("cation", cationMap, FluidComponent::Sodium, FluidComponent::MeanFieldLJ, false)
484 { require("fluid-solvent"); //which in turn requires fluid indirectly
485 }
486
processCommandFluidCation487 void process(ParamList& pl, Everything& e)
488 { CommandFluidComponent::process(pl, e);
489 }
490
printStatusCommandFluidCation491 void printStatus(Everything& e, int iRep)
492 { print(e, *(e.eVars.fluidParams.cations[iRep]));
493 }
494 }
495 commandFluidCation;
496
497 struct CommandFluidAnion : public CommandFluidComponent
498 {
CommandFluidAnionCommandFluidAnion499 CommandFluidAnion() : CommandFluidComponent("anion", anionMap, FluidComponent::Chloride, FluidComponent::MeanFieldLJ, false)
500 { require("fluid-solvent"); //which in turn requires fluid indirectly
501 }
502
processCommandFluidAnion503 void process(ParamList& pl, Everything& e)
504 { CommandFluidComponent::process(pl, e);
505 }
506
printStatusCommandFluidAnion507 void printStatus(Everything& e, int iRep)
508 { print(e, *(e.eVars.fluidParams.anions[iRep]));
509 }
510 }
511 commandFluidAnion;
512
513 enum FluidSiteParameter
514 {
515 FSp_Znuc, //!<magnitude of the nuclear charge (positive)
516 FSp_sigmaNuc, //!<gaussian width of the nuclear charge (positive)
517 FSp_Zelec, //!<magnitude of electron charge (positive)
518 FSp_aElec, //!<exponential decay width of electron charge distribution
519 FSp_sigmaElec, //!<width of peak in electron charge distribution
520 FSp_rcElec, //!<Location of peak in electron charge distribution
521 FSp_elecFilename, //!<filename to read in additional radial realspace electron charge distribution,
522 FSp_elecFilenameG, //!<filename to read in additional radial Gspace electron charge distribution
523 FSp_alpha, //!<isotropic polarizability
524 FSp_aPol, //!<cuspless-exponential width of polarizability
525 FSp_Rhs, //!< hard sphere radius
526 FSp_Delim //!< Delimiter used in parsing:
527 };
528
529 EnumStringMap<FluidSiteParameter> FSParamMap(
530 FSp_Znuc, "Znuc",
531 FSp_sigmaNuc, "sigmaNuc",
532 FSp_Zelec, "Zelec",
533 FSp_aElec, "aElec",
534 FSp_sigmaElec, "sigmaElec",
535 FSp_rcElec,"rcElec",
536 FSp_alpha, "alpha",
537 FSp_aPol, "aPol",
538 FSp_Rhs, "Rhs",
539 FSp_elecFilename, "elecFilename",
540 FSp_elecFilenameG, "elecFilenameG" );
541
542 EnumStringMap<FluidSiteParameter> FSParamDescMap(
543 FSp_Znuc, "magnitude of the nuclear charge (positive)",
544 FSp_sigmaNuc, "gaussian width of the nuclear charge (positive)",
545 FSp_Zelec, "magnitude of electron charge (positive)",
546 FSp_aElec, "exponential decay width of electron charge distribution",
547 FSp_sigmaElec, "width of peak in electron charge distribution",
548 FSp_rcElec, "location of peak in electron charge distribution",
549 FSp_elecFilename, "filename to read in additional radial realspace electron charge distribution",
550 FSp_elecFilenameG, "filename to read in additional radial Gspace electron charge distribution",
551 FSp_alpha, "isotropic polarizability",
552 FSp_aPol, "cuspless-exponential width of polarizability",
553 FSp_Rhs, "hard sphere radius for use in FMT" );
554
555 //Kendra: list of fluid components supported
556 EnumStringMap<FluidComponent::Name> fluidComponentMap(
557 //solvents
558 FluidComponent::H2O, "H2O",
559 FluidComponent::CHCl3, "CHCl3",
560 FluidComponent::CCl4, "CCl4",
561 FluidComponent::CH3CN, "CH3CN",/*
562 FluidComponent::DMC, "DMC",
563 FluidComponent::EC, "EC",
564 FluidComponent::PC, "PC",
565 FluidComponent::DMF, "DMF",
566 FluidComponent::THF, "THF",
567 FluidComponent::EthylEther, "EthylEther",
568 FluidComponent::Chlorobenzene, "Chlorobenzene",
569 FluidComponent::Isobutanol, "Isobutanol",
570 FluidComponent::CarbonDisulfide, "CarbonDisulfide",
571 FluidComponent::CustomSolvent, "CustomSolvent",*/
572 //cations
573 FluidComponent::Sodium, "Na+",
574 FluidComponent::HydratedSodium, "Na(H2O)4+",
575 FluidComponent::CustomCation, "CustomCation",
576 //anions
577 FluidComponent::Chloride, "Cl-",
578 FluidComponent::Fluoride, "F-",
579 FluidComponent::Perchlorate, "ClO4-",
580 FluidComponent::CustomAnion, "CustomAnion" );
581
582 struct CommandFluidSiteParams : public Command
583 {
CommandFluidSiteParamsCommandFluidSiteParams584 CommandFluidSiteParams() : Command("fluid-site-params", "jdftx/Fluid/Constituents")
585 {
586 format = " <component> <siteName> <key1> <value1> <key2> <value2> ...";
587 comments = "Set parameters of site <siteName> for fluid <component> which may be one of:"
588 + addDescriptions(fluidComponentMap.optionList(), nullDescription)
589 + "\n\nPossible keys and value types are:"
590 + addDescriptions(FSParamMap.optionList(), linkDescription(FSParamMap, FSParamDescMap))
591 + "\n\nAny number of these key-value pairs may be specified in any order.";
592
593 require("fluid-solvent");
594 hasDefault = false;
595 allowMultiple = true;
596 }
597
processCommandFluidSiteParams598 void process(ParamList& pl, Everything& e)
599 {
600 if(e.eVars.fluidParams.fluidType == FluidNone)
601 return;
602 FluidSolverParams& fsp = e.eVars.fluidParams;
603
604 //Read in and check name of the solvent, get index of the solvent in FluidComponent
605 FluidComponent::Name solventName;
606 pl.get(solventName, fsp.components[0]->name, fluidComponentMap, "solvent", false);
607 std::shared_ptr<FluidComponent> FC;
608 for (const auto& c : fsp.components)
609 {
610 if (solventName == c->name)
611 FC=c;
612 }
613 if (!FC)
614 throw string("Choice of <solvent> is not valid.\n Hint: Issue fluid-solvent first");
615
616 //Read in name of the site
617 string siteName;
618 pl.get(siteName, FC->molecule.sites[0]->name, "siteName", false);
619 std::shared_ptr<Molecule::Site> site;
620 for (const auto& s : FC->molecule.sites)
621 {
622 if(siteName == s->name)
623 site=s;
624 }
625 if (!site)
626 throw string("Choice of <siteName> is not valid.");
627
628 //Read parameters:
629 while(true)
630 { FluidSiteParameter key;
631 pl.get(key, FSp_Delim, FSParamMap, "key");
632 #define READ_AND_CHECK(param, op, val) \
633 case FSp_##param: \
634 pl.get(site->param, val, #param, true); \
635 if(!(site->param op val)) throw string(#param " must be " #op " " #val); \
636 break;
637 switch(key)
638 {
639 READ_AND_CHECK(Znuc,>=,0.)
640 READ_AND_CHECK(sigmaNuc,>=,0.)
641 READ_AND_CHECK(Zelec,>=,0.)
642 READ_AND_CHECK(aElec,>,0.)
643 READ_AND_CHECK(sigmaElec,>=,0.)
644 READ_AND_CHECK(rcElec,>=,0.)
645 READ_AND_CHECK(elecFilename,!=,string(""))
646 READ_AND_CHECK(elecFilenameG,!=,string(""))
647 READ_AND_CHECK(alpha,>,0.)
648 READ_AND_CHECK(aPol,>,0.)
649 READ_AND_CHECK(Rhs,>,0.)
650 case FSp_Delim: return; //end of input
651 }
652 #undef READ_AND_CHECK
653 }
654 }
655
printStatusCommandFluidSiteParams656 void printStatus(Everything& e, int iRep)
657 {
658 //prints all the sites and parameters, even if the default is unchanged
659 #define PRINT(param) logPrintf(" \\\n\t" #param " %lg", s->param);
660
661 if(e.eVars.fluidParams.fluidType == FluidNonlinearPCM || e.eVars.fluidParams.fluidType == FluidLinearPCM || e.eVars.fluidParams.fluidType == FluidNone)
662 return;
663 if(iRep==0)
664 {
665 int counter=0;
666 const FluidSolverParams& fsp = e.eVars.fluidParams;
667 for (const auto& c : fsp.components)
668 {
669 string cName = fluidComponentMap.getString(c->name);
670 for (const auto& s : c->molecule.sites)
671 {
672 string sName = s->name;
673 if(counter)
674 logPrintf("\nfluid-site-params ");
675 logPrintf("%s %s",cName.c_str(),sName.c_str()),
676 #define PRINT(param) logPrintf(" \\\n\t" #param " %lg", s->param);
677 PRINT(Znuc)
678 PRINT(sigmaNuc)
679 PRINT(Zelec)
680 PRINT(aElec)
681 PRINT(sigmaElec)
682 PRINT(rcElec)
683 PRINT(alpha)
684 PRINT(aPol)
685 PRINT(Rhs)
686 logPrintf(" \\\n\telecFilename ");
687 if (s->elecFilename.length())
688 logPrintf("%s", s->elecFilename.c_str());
689 logPrintf(" \\\n\telecFilenameG ");
690 if (s->elecFilenameG.length())
691 logPrintf("%s", s->elecFilenameG.c_str());
692 #undef PRINT
693
694 counter++;
695 }
696 }
697 }
698 }
699 }
700 commandFSParams;
701
702
703 EnumStringMap<FMixFunctional> fMixMap
704 (
705 LJPotential, "LJPotential",
706 GaussianKernel, "GaussianKernel"
707 );
708
709 struct CommandFluidMixingFunctional : public Command
710 {
CommandFluidMixingFunctionalCommandFluidMixingFunctional711 CommandFluidMixingFunctional() : Command("fluid-mixing-functional", "jdftx/Fluid/Constituents")
712 {
713 format = "<fluid1> <fluid2> <energyScale> [<lengthScale>] [<FMixType>=LJPotential]";
714 comments =
715 "Couple named fluids <fluid1> and <fluid2> which could each be one of:"
716 + addDescriptions(fluidComponentMap.optionList(), nullDescription)
717 + "\n\ntogether with a mixing functional of type:"
718 + addDescriptions(fMixMap.optionList(), nullDescription)
719 + "\n\nwith strength <energyScale> (in Eh) and range parameter <lengthScale> (in bohrs).\n";
720
721 require("fluid-solvent"); //which in turn requires fluid indirectly
722 allowMultiple = true;
723 }
724
processCommandFluidMixingFunctional725 void process(ParamList& pl, Everything& e)
726 {
727 FluidSolverParams& fsp = e.eVars.fluidParams;
728 FmixParams fmp;
729
730 FluidComponent::Name name1;
731 pl.get(name1, fsp.components[0]->name, fluidComponentMap, "fluid1", true);
732
733 FluidComponent::Name name2;
734 pl.get(name2, fsp.components[0]->name, fluidComponentMap, "fluid2", true);
735
736
737 for(const std::shared_ptr<FluidComponent> c: fsp.components)
738 {
739 if(c->name == name1) fmp.fluid1 = c;
740 if(c->name == name2) fmp.fluid2 = c;
741 }
742
743 if(!fmp.fluid1)
744 throw string("Choice of <fluid1> = %s is not valid.\n Hint: Issue fluid-solvent first.",name1);
745 if(!fmp.fluid2)
746 throw string("Choice of <fluid2> = %s is not valid.\n Hint: Issue fluid-solvent first.",name2);
747 if(fmp.fluid1->name == fmp.fluid2->name)
748 throw string("<fluid1>=<fluid2> Cannot specify mixing functional for the same fluid.");
749
750 double default_energyscale = sqrt(fmp.fluid1->epsLJ*fmp.fluid2->epsLJ);
751 pl.get(fmp.energyScale, default_energyscale,"energyScale", true);
752
753 double default_lengthscale = fmp.fluid1->Rvdw + fmp.fluid2->Rvdw;
754 pl.get(fmp.lengthScale, default_lengthscale,"lengthScale");
755
756 FMixFunctional defaultFunctional = LJPotential;
757 pl.get(fmp.FmixType, defaultFunctional, fMixMap, "FMixType");
758
759 fsp.FmixList.push_back(fmp);
760 }
761
printStatusCommandFluidMixingFunctional762 void printStatus(Everything& e, int iRep)
763 {
764 const FluidSolverParams& fsp = e.eVars.fluidParams;
765 const FmixParams& fmp = fsp.FmixList[iRep];
766
767 string c1Name = fluidComponentMap.getString(fmp.fluid1->name);
768 string c2Name = fluidComponentMap.getString(fmp.fluid2->name);
769 string fmixName = fMixMap.getString(fmp.FmixType);
770
771 logPrintf("%s %s %lg %lg %s",c1Name.c_str(),c2Name.c_str(),fmp.energyScale,fmp.lengthScale, fmixName.c_str());
772 }
773 }
774 commandFluidMixingFunctional;
775
776
777 struct CommandFluidDielectricConstant : public Command
778 {
CommandFluidDielectricConstantCommandFluidDielectricConstant779 CommandFluidDielectricConstant() : Command("fluid-dielectric-constant", "jdftx/Fluid/Parameters")
780 {
781 format = "[<epsBulkOverride>=0] [<epsInfOverride>=0]";
782 comments = "Override bulk static or high frequency dieelctric constant of fluid (if non-zero values specified)";
783 }
784
processCommandFluidDielectricConstant785 void process(ParamList& pl, Everything& e)
786 { FluidSolverParams& fsp = e.eVars.fluidParams;
787 pl.get(fsp.epsBulkOverride, 0., "epsBulkOverride");
788 pl.get(fsp.epsInfOverride, 0., "epsInfOverride");
789 }
790
printStatusCommandFluidDielectricConstant791 void printStatus(Everything& e, int iRep)
792 { const FluidSolverParams& fsp = e.eVars.fluidParams;
793 logPrintf("%lg %lg", fsp.epsBulkOverride, fsp.epsInfOverride);
794 }
795 }
796 commandFluidDielectricConstant;
797
798 struct CommandFluidDielectricTensor : public Command
799 {
CommandFluidDielectricTensorCommandFluidDielectricTensor800 CommandFluidDielectricTensor() : Command("fluid-dielectric-tensor", "jdftx/Fluid/Parameters")
801 {
802 format = "<epsBulkXX> <epsBulkYY> <epsBulkZZ>";
803 comments =
804 "Override bulk static dielectric constant of fluid with a tensor, assuming\n"
805 "that the Cartesian axes are the principal axes, without loss of generality.\n"
806 "Supported only for LinearPCM.";
807 require("fluid");
808 }
809
processCommandFluidDielectricTensor810 void process(ParamList& pl, Everything& e)
811 { FluidSolverParams& fsp = e.eVars.fluidParams;
812 pl.get(fsp.epsBulkTensor[0], 0., "epsBulkXX");
813 pl.get(fsp.epsBulkTensor[1], 0., "epsBulkYY");
814 pl.get(fsp.epsBulkTensor[2], 0., "epsBulkZZ");
815 if(fsp.fluidType != FluidLinearPCM)
816 throw string("Anisotropic epsilon supported only for LinearPCM.");
817 }
818
printStatusCommandFluidDielectricTensor819 void printStatus(Everything& e, int iRep)
820 { const vector3<>& eps = e.eVars.fluidParams.epsBulkTensor;
821 logPrintf("%lg %lg %lg", eps[0], eps[1], eps[2]);
822 }
823 }
824 commandFluidDielectricTensor;
825