1 /*-------------------------------------------------------------------
2 Copyright 2013 Ravishankar Sundararaman, Kendra Letchworth Weaver, Deniz Gunceler
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 <fluid/FluidSolverParams.h>
21 #include <electronic/SpeciesInfo.h>
22 #include <core/Units.h>
23
FluidSolverParams()24 FluidSolverParams::FluidSolverParams()
25 : T(298*Kelvin), P(1.01325*Bar), epsBulkOverride(0.), epsInfOverride(0.), verboseLog(false), solveFrequency(FluidFreqDefault),
26 components(components_), solvents(solvents_), cations(cations_), anions(anions_),
27 vdwScale(0.75), pCavity(0.), lMax(3), cavityScale(1.), ionSpacing(0.),
28 zMask0(0.), zMaskH(0.), zMaskIonH(0.), zMaskSigma(0.5),
29 linearDielectric(false), linearScreening(false), nonlinearSCF(false), screenOverride(0.)
30 {
31 }
32
addComponent(const std::shared_ptr<FluidComponent> & component)33 void FluidSolverParams::addComponent(const std::shared_ptr<FluidComponent>& component)
34 { components_.push_back(component);
35 switch(component->type)
36 { case FluidComponent::Solvent: solvents_.push_back(component); break;
37 case FluidComponent::Cation: cations_.push_back(component); break;
38 case FluidComponent::Anion: anions_.push_back(component); break;
39 }
40 }
41
42
setCDFTparams()43 void FluidSolverParams::setCDFTparams()
44 {
45 if(solvents.size() == 1)
46 {
47 switch(solvents[0]->name)
48 {
49 case FluidComponent::H2O:
50 vdwScale = 0.540;
51 break;
52 case FluidComponent::CHCl3:
53 vdwScale = 0.393;
54 break;
55 case FluidComponent::CCl4:
56 vdwScale = 0.407;
57 break;
58 case FluidComponent::CH3CN:
59 default:
60 vdwScale = 0.488;
61 initWarnings += "WARNING: Classical DFT has not been parameterized for this solvent, using 0.488 as the Van der Waals scale factor!\n";
62 break;
63
64 }
65 }
66 else
67 {
68 vdwScale = 0.488;
69 initWarnings += "WARNING: Classical DFT has not been parameterized for solvent mixtures, using 0.488 as the Van der Waals scale factor!\n";
70 }
71 }
72
setPCMparams()73 void FluidSolverParams::setPCMparams()
74 {
75 assert(solvents.size()==1);
76 const double dyn_per_cm = (1e-5*Newton)/(1e-2*meter);
77 const double GPa = 1e9*Pascal;
78
79 //Set common default values:
80 cavityTension = 0.;
81 cavityPressure = 0.;
82 cavityScale = 1.;
83 vdwScale = 1.;
84
85 //Set PCM fit parameters:
86 switch(pcmVariant)
87 { case PCM_SaLSA:
88 { nc = 1.42e-3;
89 sigma = sqrt(0.5);
90 switch(solvents[0]->name)
91 {
92 case FluidComponent::H2O:
93 vdwScale = 0.50;
94 break;
95 case FluidComponent::CHCl3:
96 vdwScale = 0.88;
97 break;
98 case FluidComponent::CCl4:
99 vdwScale = 1.06;
100 break;
101 case FluidComponent::CH3CN:
102 vdwScale = 0.37; //Tentative: using toluene solvation energy alone
103 break;
104 default:
105 vdwScale = 1.;
106 initWarnings += "WARNING: SALSA has not been parametrized for this solvent, using 1.0 as the Van der Waals scale factor!\n";
107 break;
108 }
109 assert(fluidType == FluidSaLSA);
110 break;
111 }
112 case PCM_CANDLE:
113 { nc = 1.42e-3;
114 sigma = sqrt(0.5);
115 switch(solvents[0]->name)
116 { case FluidComponent::CH3CN:
117 Ztot = 16;
118 eta_wDiel = 3.15;
119 sqrtC6eff = 2.21;
120 pCavity = -31.;
121 break;
122 case FluidComponent::Octanol:
123 Ztot = 56;
124 eta_wDiel = 5.5;
125 sqrtC6eff = 16.;
126 pCavity = 0;
127 initWarnings += "WARNING: CANDLE parameters for Octanol have not yet been fit (initial guess only)\n";
128 break;
129 case FluidComponent::DMSO:
130 Ztot = 26;
131 eta_wDiel = 3.8;
132 sqrtC6eff = 8.;
133 pCavity = 20.;
134 initWarnings += "WARNING: CANDLE parameters for DMSO have not yet been fit (initial guess only)\n";
135 break;
136 case FluidComponent::H2O:
137 default:
138 Ztot = 8;
139 eta_wDiel = 1.46;
140 sqrtC6eff = 0.770;
141 pCavity = 36.5;
142 if(solvents[0]->name != FluidComponent::H2O)
143 initWarnings += "WARNING: CANDLE LinearPCM has not been parametrized for this solvent, using fit parameters for water\n";
144 break;
145 }
146 if(fluidType != FluidLinearPCM) initWarnings += "WARNING: CANDLE has only been parametrized for LinearPCM.\n";
147 break;
148 }
149 case PCM_SGA13:
150 { nc = 1e-2;
151 sigma = 0.6;
152 switch(solvents[0]->name)
153 {
154 case FluidComponent::H2O:
155 vdwScale = 0.538;
156 break;
157 case FluidComponent::CHCl3:
158 vdwScale = 1.315;
159 break;
160 case FluidComponent::CCl4:
161 vdwScale = 1.238;
162 break;
163 case FluidComponent::CH3CN:
164 vdwScale = 0.6; //Tentative: using toluene and self solvation energies alone
165 break;
166 default:
167 vdwScale = 1.;
168 initWarnings += "WARNING: SGA13 has not been parametrized for this solvent, using 1.0 as the Van der Waals scale factor!\n";
169 break;
170 }
171 break;
172 }
173 case PCM_GLSSA13:
174 {
175 switch(solvents[0]->name)
176 {
177 case FluidComponent::CHCl3:
178 { nc = 2.4e-05;
179 sigma = 0.6;
180 cavityTension = -9.066e-6;
181 break;
182 }
183 case FluidComponent::CCl4:
184 { switch(fluidType)
185 { case FluidLinearPCM:
186 nc = 1.15e-4;
187 sigma = 0.6;
188 cavityTension = -8.99e-06;
189 break;
190 case FluidNonlinearPCM:
191 die("\nERROR: You can't use NonlinearPCM with CCl4 as it does not have a permanent dipole moment!\n");
192 default: //Other fluids do not use these parameters
193 break;
194 }
195 break;
196 }
197 case FluidComponent::CH3CN:
198 { nc = 1.8e-4;
199 sigma = 0.6;
200 cavityTension = -6.3e-7;
201 break;
202 }
203 case FluidComponent::CH2Cl2:
204 { nc = 9.3e-4;
205 sigma = 0.6;
206 cavityTension = -2.7e-6;
207 break;
208 }
209 case FluidComponent::Ethanol:
210 { nc = 1.3e-3;
211 sigma = 0.6;
212 cavityTension = -5.1e-6;
213 break;
214 }
215 case FluidComponent::Methanol:
216 { nc = 6.5e-4;
217 sigma = 0.6;
218 cavityTension = -5.2e-6;
219 break;
220 }
221 case FluidComponent::EthylEther:
222 { nc = 2.63e-4;
223 sigma = 0.6;
224 cavityTension = -1.08e-5;
225 break;
226 }
227 case FluidComponent::Chlorobenzene:
228 { nc = 4.28e-5;
229 sigma = 0.6;
230 cavityTension = -6.2e-06;
231 break;
232 }
233 case FluidComponent::Isobutanol:
234 { nc = 1.51e-3;
235 sigma = 0.6;
236 cavityTension = -8.96e-06;
237 break;
238 }
239 case FluidComponent::EC:
240 { nc = 1.8e-3;
241 sigma = 0.6;
242 cavityTension = 1.55e-5;
243 break;
244 }
245 case FluidComponent::PC:
246 { nc = 9.8e-4;
247 sigma = 0.6;
248 cavityTension = 9.53e-6;
249 break;
250 }
251 case FluidComponent::EthyleneGlycol:
252 { switch(fluidType)
253 { case FluidLinearPCM:
254 { nc = 5.4e-4;
255 sigma = 0.6;
256 cavityTension = 1.15e-5;
257 break;
258 }
259 case FluidNonlinearPCM:
260 die("\nERROR: You can't use NonlinearPCM with Ethylene Glycol as it does not have a permanent dipole moment!\n");
261 default: //Other fluids do not use these parameters
262 break;
263 }
264 }
265 case FluidComponent::Glyme:
266 { switch(fluidType)
267 { case FluidLinearPCM:
268 { nc = 8.36e-5;
269 sigma = 0.6;
270 cavityTension = -8.03e-06;
271 break;
272 }
273 case FluidNonlinearPCM:
274 die("\nERROR: You can't use NonlinearPCM with Glyme as it does not have a permanent dipole moment!\n");
275 default: //Other fluids do not use these parameters
276 break;
277 }
278 break;
279 }
280 case FluidComponent::CarbonDisulfide:
281 { switch(fluidType)
282 { case FluidLinearPCM:
283 nc = 1.6e-5;
284 sigma = 0.6;
285 cavityTension = -6.8e-06;
286 break;
287 case FluidNonlinearPCM:
288 die("\nERROR: You can't use NonlinearPCM with CarbonDisulfide as it does not have a permanent dipole moment!\n");
289 default: //Other fluids do not use these parameters
290 break;
291 }
292 break;
293 }
294 case FluidComponent::THF:
295 { nc = 1.6e-3;
296 sigma = 0.6;
297 cavityTension = -1.7e-06;
298 break;
299 }
300 case FluidComponent::DMSO:
301 { nc = 9.5e-4;
302 sigma = 0.6;
303 cavityTension = 8.42e-06;
304 break;
305 }
306 default: // For water and unparametrized fluids
307 { switch(fluidType)
308 { case FluidLinearPCM:
309 nc = 3.7e-4;
310 sigma = 0.6;
311 cavityTension = 5.4e-6;
312 break;
313 case FluidNonlinearPCM:
314 nc = 1.0e-3;
315 sigma = 0.6;
316 cavityTension = 9.5e-6;
317 break;
318 default: //Other fluids do not use these parameters
319 break;
320 }
321 if(solvents[0]->name != FluidComponent::H2O)
322 { initWarnings +=
323 "WARNING: PCM variant GLSSA has not been parametrized for this solvent; using bulk\n"
324 " surface tension as effective cavity tension and water parameters for cavity size.\n";
325 cavityTension = solvents[0]->sigmaBulk;
326 }
327 break;
328 }
329 }
330 break;
331 }
332 case PCM_LA12:
333 { nc = 7e-4;
334 sigma = 0.6;
335 if(fluidType == FluidNonlinearPCM)
336 initWarnings += "WARNING: PCM variant LA12/PRA05 has not been parametrized for NonlinearPCM; using LinearPCM fit parameters.\n";
337 if( (fluidType==FluidLinearPCM || fluidType==FluidNonlinearPCM) && solvents[0]->name != FluidComponent::H2O)
338 initWarnings += "WARNING: PCM variant LA12/PRA05 has been fit only for H2O; using nc and sigma from H2O fit.\n";
339 break;
340 }
341 case PCM_SoftSphere:
342 { sigma = 0.5; //in bohrs for this model
343 switch(solvents[0]->name)
344 { case FluidComponent::Ethanol:
345 { cavityTension = -4.0*dyn_per_cm;
346 cavityScale = 1.22;
347 break;
348 }
349 default: // For water and unparametrized fluids
350 { if(fluidType==FluidLinearPCM)
351 { cavityTension = 9.67e-6;
352 cavityScale = 1.11;
353 }
354 else
355 { cavityTension = 1.02e-5;
356 cavityScale = 1.00;
357 }
358 //Warn when not explicitly parametrized:
359 if(solvents[0]->name != FluidComponent::H2O)
360 { initWarnings +=
361 "WARNING: PCM variant SoftSphere has not been parametrized for this solvent; using bulk\n"
362 " surface tension as effective cavity tension and water parameters for everything else.\n";
363 cavityTension = solvents[0]->sigmaBulk;
364 }
365 }
366 }
367 }
368 case PCM_FixedCavity:
369 { //no default parameters
370 }
371 case_PCM_SCCS_any:
372 { if(fluidType != FluidLinearPCM) initWarnings += "WARNING: SCCS has only been parametrized for LinearPCM.\n";
373 if(solvents[0]->name != FluidComponent::H2O)
374 { initWarnings +=
375 "WARNING: SCCS variants have not been parametrized for this solvent; using water parameters\n";
376 }
377 rhoDelta = 1e-4; //common to all variants
378 switch(pcmVariant)
379 { case PCM_SCCS_g09: rhoMin=1.00e-4; rhoMax=1.50e-3; cavityTension=2.50*dyn_per_cm; break;
380 case PCM_SCCS_g03: rhoMin=1.00e-4; rhoMax=5.00e-3; cavityTension=11.5*dyn_per_cm; break;
381 case PCM_SCCS_g03p: rhoMin=3.00e-4; rhoMax=3.00e-3; cavityTension=12.0*dyn_per_cm; break;
382 case PCM_SCCS_g09beta: rhoMin=1.00e-4; rhoMax=1.50e-3; cavityTension=11.0*dyn_per_cm; cavityPressure=-0.08*GPa; break;
383 case PCM_SCCS_g03beta: rhoMin=1.00e-4; rhoMax=5.00e-3; cavityTension=50.0*dyn_per_cm; cavityPressure=-0.35*GPa; break;
384 case PCM_SCCS_g03pbeta: rhoMin=3.00e-4; rhoMax=3.00e-3; cavityTension=20.0*dyn_per_cm; cavityPressure=-0.08*GPa; break;
385 case PCM_SCCS_cation: rhoMin=2.00e-4; rhoMax=3.50e-3; cavityTension=5.00*dyn_per_cm; cavityPressure=+0.125*GPa; break;
386 case PCM_SCCS_anion: rhoMin=2.40e-3; rhoMax=1.55e-2; cavityTension=0.00*dyn_per_cm; cavityPressure=+0.450*GPa; break;
387 default: break;
388 }
389 break;
390 }
391 }
392 }
393
getAtomicRadius(const SpeciesInfo & sp) const394 double FluidSolverParams::getAtomicRadius(const SpeciesInfo& sp) const
395 { static const int Zmax = 103;
396 static double uffDiaAngst[Zmax] = {
397 2.886, 2.362, 2.451, 2.745, 4.083, 3.851, 3.100, 3.500, 3.364, 3.243, //Z -> 10 (Z=7 i.e. N overridden to Bondi vdW radius)
398 2.983, 3.021, 4.499, 4.295, 4.147, 4.035, 3.947, 3.868, 3.812, 3.399, //Z -> 20
399 3.295, 3.175, 3.144, 3.023, 2.961, 2.912, 2.872, 2.834, 3.495, 2.763, //Z -> 30
400 4.383, 4.280, 4.230, 4.205, 4.189, 4.141, 4.114, 3.641, 3.345, 3.124, //Z -> 40
401 3.165, 3.052, 2.998, 2.963, 2.929, 2.899, 3.148, 2.848, 4.463, 4.392, //Z -> 50
402 4.420, 4.470, 4.500, 4.404, 4.517, 3.703, 3.522, 3.556, 3.606, 3.575, //Z -> 60
403 3.547, 3.520, 3.493, 3.368, 3.451, 3.428, 3.409, 3.391, 3.374, 3.355, //Z -> 70
404 3.640, 3.141, 3.170, 3.069, 2.954, 3.120, 2.840, 2.754, 3.293, 2.705, //Z -> 80
405 4.347, 4.297, 4.370, 4.709, 4.750, 4.765, 4.900, 3.677, 3.478, 3.396, //Z -> 90
406 3.424, 3.395, 3.424, 3.424, 3.381, 3.326, 3.339, 3.313, 3.299, 3.286, //Z -> 100
407 3.274, 3.248, 3.236 }; //Z -> 103
408 if(sp.atomicRadiusOverride) return *(sp.atomicRadiusOverride); //Return override value, if available
409 assert(sp.atomicNumber && sp.atomicNumber <= Zmax);
410 return uffDiaAngst[sp.atomicNumber-1]*0.5*Angstrom; //convert to radius in bohrs
411 }
412
413
needsVDW() const414 bool FluidSolverParams::needsVDW() const
415 { switch(fluidType)
416 { case FluidNone:
417 return false;
418 case FluidLinearPCM:
419 case FluidNonlinearPCM:
420 return (pcmVariant==PCM_SGA13 || pcmVariant==PCM_CANDLE);
421 case FluidSaLSA:
422 case FluidClassicalDFT:
423 default:
424 return true;
425 }
426 }
427
ionicScreening() const428 bool FluidSolverParams::ionicScreening() const
429 { return cations.size() && anions.size();
430 }
431