1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman, Deniz Gunceler, 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 <commands/ParamList.h>
22 #include <electronic/Everything.h>
23 #include <electronic/IonicDynamicsParams.h>
24 #include <core/Units.h>
25
26
27 EnumStringMap<IonicDynamicsParams::StatMethod> statMethodMap
28 ( IonicDynamicsParams::StatNone, "None",
29 IonicDynamicsParams::Berendsen, "Berendsen",
30 IonicDynamicsParams::NoseHoover, "NoseHoover"
31 );
32
33 //An enum entry for each configurable option of IonicDynamicsParams
34 enum IonicDynamicsParamsMember
35 { IDPM_dt,
36 IDPM_nSteps,
37 IDPM_statMethod,
38 IDPM_T0,
39 IDPM_P0,
40 IDPM_stress0,
41 IDPM_tDampT,
42 IDPM_tDampP,
43 IDPM_chainLengthT,
44 IDPM_chainLengthP,
45 IDPM_B0,
46 IDPM_Delim //!< delimiter to detect end of input
47 };
48
49 EnumStringMap<IonicDynamicsParamsMember> idpmMap
50 ( IDPM_dt, "dt",
51 IDPM_nSteps, "nSteps",
52 IDPM_statMethod, "statMethod",
53 IDPM_T0, "T0",
54 IDPM_P0, "P0",
55 IDPM_stress0, "stress0",
56 IDPM_tDampT, "tDampT",
57 IDPM_tDampP, "tDampP",
58 IDPM_chainLengthT, "chainLengthT",
59 IDPM_chainLengthP, "chainLengthP",
60 IDPM_B0, "B0"
61 );
62
63 EnumStringMap<IonicDynamicsParamsMember> idpmDescMap
64 ( IDPM_dt, "time step [fs]",
65 IDPM_nSteps, "number of molecular dynamics steps",
66 IDPM_statMethod, statMethodMap.optionList() + " (method for thermostat and/or barostat)",
67 IDPM_T0, "thermostat temperature or temperature for initial velocities [Kelvin]",
68 IDPM_P0, "barostat pressure [bar] (default NAN => no hydrostatic barostat)",
69 IDPM_stress0, "barostat stress components xx yy zz yz zx xy [bar] (default NANs => no anisotropic barostat)",
70 IDPM_tDampT, "thermostat damping time [fs]",
71 IDPM_tDampP, "barostat damping time [fs]",
72 IDPM_chainLengthT, "Nose-Hoover chain length for thermostat",
73 IDPM_chainLengthP, "Nose-Hoover chain length for barostat",
74 IDPM_B0, "Characteristic bulk modulus [bar] for Berendsen barostat (damping ~ B0 * tDampP)"
75 );
76
77 struct CommandIonicDynamics : public Command
78 {
CommandIonicDynamicsCommandIonicDynamics79 CommandIonicDynamics() : Command("ionic-dynamics", "jdftx/Ionic/Optimization")
80 { format = "<time-step> <total-time> [<kT>=0.001] [<alpha>=0.0]";
81 comments = "Applies Verlet molecular dynamics\n"
82 "Takes the time in femtoseconds and kT is in Hartree.\n"
83 "Give <alpha> if you want to the system to equilibrate with a heat bath at <kT>.\n"
84 "Also make sure to turn the symmetries off if the initial velocities don't satisfy the symmetry\n"
85 "conditions. Initial velocities will be assigned randomly if they are not given.";
86
87 format = "<key1> <value1> <key2> <value2> ...";
88 comments = "Born-Oppenheimer molecular dynamics, controlled by keys:"
89 + addDescriptions(idpmMap.optionList(), linkDescription(idpmMap, idpmDescMap))
90 + "\n\nAny number of these key-value pairs may be specified in any order.\n\n"
91 "Note that nSteps must be non-zero to activate dynamics.\n"
92 "Default mode is NVE; specify statMethod to add a thermostat or barostat";
93 }
94
processCommandIonicDynamics95 void process(ParamList& pl, Everything& e)
96 { IonicDynamicsParams& idp = e.ionicDynParams;
97 const double nanVal = NAN;
98 while(true)
99 { IonicDynamicsParamsMember key;
100 pl.get(key, IDPM_Delim, idpmMap, "key");
101 switch(key)
102 { case IDPM_dt: pl.get(idp.dt, 1., "dt", true); idp.dt *= fs; break;
103 case IDPM_nSteps: pl.get(idp.nSteps, 0, "nSteps", true); break;
104 case IDPM_statMethod: pl.get(idp.statMethod, IonicDynamicsParams::StatNone, statMethodMap, "statMethod", true); break;
105 case IDPM_T0: pl.get(idp.T0, nanVal, "T0", true); idp.T0 *= Kelvin; break;
106 case IDPM_P0: pl.get(idp.P0, nanVal, "P0", true); idp.P0 *= Bar; break;
107 case IDPM_stress0:
108 { pl.get(idp.stress0(0,0), nanVal, "stress0_xx", true);
109 pl.get(idp.stress0(1,1), nanVal, "stress0_yy", true);
110 pl.get(idp.stress0(2,2), nanVal, "stress0_zz", true);
111 pl.get(idp.stress0(1,2), nanVal, "stress0_yz", true);
112 pl.get(idp.stress0(2,0), nanVal, "stress0_zx", true);
113 pl.get(idp.stress0(0,1), nanVal, "stress0_xy", true);
114 idp.stress0(2,1) = idp.stress0(1,2);
115 idp.stress0(0,2) = idp.stress0(2,0);
116 idp.stress0(1,0) = idp.stress0(0,1);
117 idp.stress0 *= Bar;
118 break;
119 }
120 case IDPM_tDampT: pl.get(idp.tDampT, 50., "tDampT", true); idp.tDampT *= fs; break;
121 case IDPM_tDampP: pl.get(idp.tDampP, 100., "tDampP", true); idp.tDampP *= fs; break;
122 case IDPM_chainLengthT: pl.get(idp.chainLengthT, 3, "chainLengthT", true); break;
123 case IDPM_chainLengthP: pl.get(idp.chainLengthP, 3, "chainLengthP", true); break;
124 case IDPM_B0: pl.get(idp.B0, nanVal, "B0", true); idp.B0 *= Bar; break;
125 case IDPM_Delim:
126 if((not std::isnan(idp.P0)) and (not std::isnan(trace(idp.stress0))))
127 throw(string("Cannot specify both P0 (hydrostatic) and stress0 (anisotropic) barostats"));
128 return; //end of input
129 }
130 }
131 }
132
printStatusCommandIonicDynamics133 void printStatus(Everything& e, int iRep)
134 { const IonicDynamicsParams& idp = e.ionicDynParams;
135 logPrintf(" \\\n\tdt %lg", idp.dt/fs);
136 logPrintf(" \\\n\tnSteps %d", idp.nSteps);
137 logPrintf(" \\\n\tstatMethod %s", statMethodMap.getString(idp.statMethod));
138 logPrintf(" \\\n\tT0 %lg", idp.T0/Kelvin);
139 logPrintf(" \\\n\tP0 %lg", idp.P0/Bar);
140 logPrintf(" \\\n\tstress0 %lg %lg %lg %lg %lg %lg",
141 idp.stress0(0,0)/Bar, idp.stress0(1,1)/Bar, idp.stress0(2,2)/Bar,
142 idp.stress0(1,2)/Bar, idp.stress0(2,0)/Bar, idp.stress0(0,1)/Bar);
143 logPrintf(" \\\n\ttDampT %lg", idp.tDampT/fs);
144 logPrintf(" \\\n\ttDampP %lg", idp.tDampP/fs);
145 logPrintf(" \\\n\tchainLengthT %d", idp.chainLengthT);
146 logPrintf(" \\\n\tchainLengthP %d", idp.chainLengthP);
147 logPrintf(" \\\n\tB0 %lg", idp.B0/Bar);
148 }
149 }
150 commandIonicDynamics;
151
152
153 struct CommandLjOverride : public Command
154 {
CommandLjOverrideCommandLjOverride155 CommandLjOverride() : Command("lj-override", "jdftx/Ionic/Optimization")
156 { format = "<rCut>";
157 comments =
158 "Replace electronic DFT by a Lennard-Jones only potential with cutoff <rCut> in Angstroms.\n"
159 "This potential will use LJ parameters, sigma = 2^(5/6) R0 and epsilon = C6/(128 R0^6),\n"
160 "where R0 and C6 are DFT-D2 parameters for each species (adjustable by command setVDW).\n"
161 "This is not for real calculations, but a quick way to debug ionic-minimize,\n"
162 "lattice-minimize or ionic-dynamics. Tip: set elec-cutoff to a small value to\n"
163 "speed up electronic initialization (which is not bypassed for code simplicity).";
164 }
165
processCommandLjOverride166 void process(ParamList& pl, Everything& e)
167 { pl.get(e.iInfo.ljOverride, 0., "rCut", true);
168 e.iInfo.ljOverride *= Angstrom; //convert to atomic units
169 }
170
printStatusCommandLjOverride171 void printStatus(Everything& e, int iRep)
172 { logPrintf("%lg", e.iInfo.ljOverride/Angstrom);
173 }
174 }
175 commandLjOverride;
176
177 //---- Thermostat / barostat velocities ----
178 struct CommandStatVelocity : public Command
179 {
CommandStatVelocityCommandStatVelocity180 CommandStatVelocity(string statName) : Command(statName+"-velocity", "jdftx/Ionic/Optimization")
181 { format = "<v1> <v2> ...";
182 comments =
183 "Read "+statName+" internal velocities for continuing ionic dynamics.\n"
184 "This command is automatically dumped with ionpos from dynamics simulations\n"
185 "using Nose-Hoover chains that involve "+statName+" internal velocities.";
186 }
187
188 virtual diagMatrix& target(Everything& e)=0;
189
processCommandStatVelocity190 void process(ParamList& pl, Everything& e)
191 { diagMatrix& stat = target(e);
192 stat.clear();
193 while(true)
194 { double v = NAN;
195 pl.get(v, v, "v");
196 if(std::isnan(v)) break;
197 stat.push_back(v);
198 }
199 }
200
printStatusCommandStatVelocity201 void printStatus(Everything& e, int iRep)
202 { const diagMatrix& stat = target(e);
203 for(const double& v: stat) logPrintf("%lg ", v);
204 }
205 };
206
207 struct CommandThermostatVelocity : public CommandStatVelocity
CommandThermostatVelocityCommandThermostatVelocity208 { CommandThermostatVelocity() : CommandStatVelocity("thermostat") {}
209 diagMatrix& target(Everything& e) { return e.iInfo.thermostat; }
210 }
211 commandThermostatVelocity;
212
213 struct CommandBarostatVelocity : public CommandStatVelocity
CommandBarostatVelocityCommandBarostatVelocity214 { CommandBarostatVelocity() : CommandStatVelocity("barostat")
215 { comments += "\n(The first six components are strain rate, while the rest are lattice thermostat velocities.)\n";
216 }
217 diagMatrix& target(Everything& e) { return e.iInfo.barostat; }
218 }
219 commandBarostatVelocity;
220