source: src/config.cpp@ 62f793

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since 62f793 was 62f793, checked in by Frederik Heber <heber@…>, 17 years ago

thermostats enumerator, necessary variables included in class config, molecule::Thermostats() and ..::Constrained
Potential()

Thermostat header definitions implemented. Can be parsed from the ESPACK config file into class config
ConstrainedPotential() calculates the transformation between two given configurations my minimsing a penalty func
tional of the distance to travel per atom. This was implemented due to Michael Griebel's talk during the MultiMat

Closing meeting in order to produce some visuals. It basically mimics a "Bain Transformation". But it is not yet
perfectly implemented.

Also, MD was corrected in VerletIntegration(). However, forces are still wrong with mpqc, as the kinetic energy increases dramatically during the MD simulation.
  • Property mode set to 100644
File size: 70.4 KB
Line 
1/** \file config.cpp
2 *
3 * Function implementations for the class config.
4 *
5 */
6
7#include "molecules.hpp"
8
9/************************************* Functions for class config ***************************/
10
11/** Constructor for config file class.
12 */
13config::config()
14{
15 mainname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
16 defaultpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
17 pseudopotpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
18 configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
19 configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
20 ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
21 ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
22 for (int j=0;j<MaxThermostats;j++)
23 ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
24 Thermostat = 4;
25 alpha = 0.;
26 ScaleTempStep = 25;
27 TempFrequency = 2.5;
28 strcpy(mainname,"pcp");
29 strcpy(defaultpath,"not specified");
30 strcpy(pseudopotpath,"not specified");
31 configpath[0]='\0';
32 configname[0]='\0';
33
34 strcpy(ThermostatNames[0],"None");
35 ThermostatImplemented[0] = 1;
36 strcpy(ThermostatNames[1],"Woodcock");
37 ThermostatImplemented[1] = 1;
38 strcpy(ThermostatNames[2],"Gaussian");
39 ThermostatImplemented[2] = 1;
40 strcpy(ThermostatNames[3],"Langevin");
41 ThermostatImplemented[3] = 1;
42 strcpy(ThermostatNames[4],"Berendsen");
43 ThermostatImplemented[4] = 1;
44 strcpy(ThermostatNames[5],"NoseHoover");
45 ThermostatImplemented[5] = 1;
46
47 FastParsing = false;
48 ProcPEGamma=8;
49 ProcPEPsi=1;
50 DoOutVis=0;
51 DoOutMes=1;
52 DoOutNICS=0;
53 DoOutOrbitals=0;
54 DoOutCurrent=0;
55 DoPerturbation=0;
56 DoFullCurrent=0;
57 DoWannier=0;
58 DoConstrainedMD=0;
59 CommonWannier=0;
60 SawtoothStart=0.01;
61 VectorPlane=0;
62 VectorCut=0;
63 UseAddGramSch=1;
64 Seed=1;
65
66 MaxOuterStep=0;
67 Deltat=1;
68 OutVisStep=10;
69 OutSrcStep=5;
70 TargetTemp=0.00095004455;
71 ScaleTempStep=25;
72 MaxPsiStep=0;
73 EpsWannier=1e-7;
74
75 MaxMinStep=100;
76 RelEpsTotalEnergy=1e-7;
77 RelEpsKineticEnergy=1e-5;
78 MaxMinStopStep=1;
79 MaxMinGapStopStep=0;
80 MaxInitMinStep=100;
81 InitRelEpsTotalEnergy=1e-5;
82 InitRelEpsKineticEnergy=1e-4;
83 InitMaxMinStopStep=1;
84 InitMaxMinGapStopStep=0;
85
86 //BoxLength[NDIM*NDIM];
87
88 ECut=128.;
89 MaxLevel=5;
90 RiemannTensor=0;
91 LevRFactor=0;
92 RiemannLevel=0;
93 Lev0Factor=2;
94 RTActualUse=0;
95 PsiType=0;
96 MaxPsiDouble=0;
97 PsiMaxNoUp=0;
98 PsiMaxNoDown=0;
99 AddPsis=0;
100
101 RCut=20.;
102 StructOpt=0;
103 IsAngstroem=1;
104 RelativeCoord=0;
105 MaxTypes=0;
106};
107
108
109/** Destructor for config file class.
110 */
111config::~config()
112{
113 Free((void **)&mainname, "config::~config: *mainname");
114 Free((void **)&defaultpath, "config::~config: *defaultpath");
115 Free((void **)&pseudopotpath, "config::~config: *pseudopotpath");
116 Free((void **)&configpath, "config::~config: *configpath");
117 Free((void **)&configname, "config::~config: *configname");
118};
119
120/** Readin of Thermostat related values from parameter file.
121 * \param *source parameter file
122 */
123void config::InitThermostats(ifstream *source)
124{
125 char *thermo = MallocString(12, "IonsInitRead: thermo");
126 int verbose = 0;
127
128 // read desired Thermostat from file along with needed additional parameters
129 if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
130 if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
131 if (ThermostatImplemented[0] == 1) {
132 Thermostat = None;
133 } else {
134 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
135 Thermostat = None;
136 }
137 } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
138 if (ThermostatImplemented[1] == 1) {
139 Thermostat = Woodcock;
140 ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
141 } else {
142 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
143 Thermostat = None;
144 }
145 } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
146 if (ThermostatImplemented[2] == 1) {
147 Thermostat = Gaussian;
148 ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
149 } else {
150 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
151 Thermostat = None;
152 }
153 } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
154 if (ThermostatImplemented[3] == 1) {
155 Thermostat = Langevin;
156 ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
157 if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
158 cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
159 } else {
160 alpha = 1.;
161 }
162 } else {
163 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
164 Thermostat = None;
165 }
166 } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
167 if (ThermostatImplemented[4] == 1) {
168 Thermostat = Berendsen;
169 ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
170 } else {
171 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
172 Thermostat = None;
173 }
174 } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
175 if (ThermostatImplemented[5] == 1) {
176 Thermostat = NoseHoover;
177 ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
178 alpha = 0.;
179 } else {
180 cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
181 Thermostat = None;
182 }
183 } else {
184 cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
185 Thermostat = None;
186 }
187 } else {
188 if ((MaxOuterStep > 0) && (TargetTemp != 0))
189 cout << Verbose(2) << "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
190 Thermostat = None;
191 }
192 Free((void **)&thermo, "InitThermostats: thermo");
193};
194
195
196/** Displays menu for editing each entry of the config file.
197 * Nothing fancy here, just lots of cout << Verbose(0)s for the menu and a switch/case
198 * for each entry of the config file structure.
199 */
200void config::Edit(molecule *mol)
201{
202 char choice;
203
204 do {
205 cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
206 cout << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl;
207 cout << Verbose(0) << " B - Default path (for runtime files)" << endl;
208 cout << Verbose(0) << " C - Path of pseudopotential files" << endl;
209 cout << Verbose(0) << " D - Number of coefficient sharing processes" << endl;
210 cout << Verbose(0) << " E - Number of wave function sharing processes" << endl;
211 cout << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl;
212 cout << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl;
213 cout << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl;
214 cout << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl;
215 cout << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl;
216 cout << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl;
217 cout << Verbose(0) << " L - 0: Wannier centres as calculated, 1: common centre for all, 2: unite centres according to spread, 3: cell centre, 4: shifted to nearest grid point" << endl;
218 cout << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl;
219 cout << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl;
220 cout << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl;
221 cout << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl;
222 cout << Verbose(0) << " Q - Initial integer value of random number generator" << endl;
223 cout << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl;
224 cout << Verbose(0) << " T - Output visual after ...th step" << endl;
225 cout << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl;
226 cout << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl;
227 cout << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl;
228 cout << Verbose(0) << " Z - Maximum number of minimization iterations" << endl;
229 cout << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl;
230 cout << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl;
231 cout << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl;
232 cout << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl;
233 cout << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl;
234 cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
235 cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
236 cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
237 cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
238 cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
239 cout << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl;
240 cout << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl;
241 cout << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl;
242 cout << Verbose(0) << " p - Number of Riemann levels" << endl;
243 cout << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl;
244 cout << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl;
245 cout << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl;
246 cout << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl;
247 cout << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl;
248 cout << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl;
249 cout << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl;
250 cout << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl;
251 cout << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl;
252 cout << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl;
253 cout << Verbose(0) << "=========================================================" << endl;
254 cout << Verbose(0) << "INPUT: ";
255 cin >> choice;
256
257 switch (choice) {
258 case 'A': // mainname
259 cout << Verbose(0) << "Old: " << config::mainname << "\t new: ";
260 cin >> config::mainname;
261 break;
262 case 'B': // defaultpath
263 cout << Verbose(0) << "Old: " << config::defaultpath << "\t new: ";
264 cin >> config::defaultpath;
265 break;
266 case 'C': // pseudopotpath
267 cout << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ";
268 cin >> config::pseudopotpath;
269 break;
270
271 case 'D': // ProcPEGamma
272 cout << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ";
273 cin >> config::ProcPEGamma;
274 break;
275 case 'E': // ProcPEPsi
276 cout << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ";
277 cin >> config::ProcPEPsi;
278 break;
279 case 'F': // DoOutVis
280 cout << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ";
281 cin >> config::DoOutVis;
282 break;
283 case 'G': // DoOutMes
284 cout << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ";
285 cin >> config::DoOutMes;
286 break;
287 case 'H': // DoOutOrbitals
288 cout << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ";
289 cin >> config::DoOutOrbitals;
290 break;
291 case 'I': // DoOutCurrent
292 cout << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ";
293 cin >> config::DoOutCurrent;
294 break;
295 case 'J': // DoFullCurrent
296 cout << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ";
297 cin >> config::DoFullCurrent;
298 break;
299 case 'K': // DoPerturbation
300 cout << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ";
301 cin >> config::DoPerturbation;
302 break;
303 case 'L': // CommonWannier
304 cout << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ";
305 cin >> config::CommonWannier;
306 break;
307 case 'M': // SawtoothStart
308 cout << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ";
309 cin >> config::SawtoothStart;
310 break;
311 case 'N': // VectorPlane
312 cout << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ";
313 cin >> config::VectorPlane;
314 break;
315 case 'O': // VectorCut
316 cout << Verbose(0) << "Old: " << config::VectorCut << "\t new: ";
317 cin >> config::VectorCut;
318 break;
319 case 'P': // UseAddGramSch
320 cout << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ";
321 cin >> config::UseAddGramSch;
322 break;
323 case 'Q': // Seed
324 cout << Verbose(0) << "Old: " << config::Seed << "\t new: ";
325 cin >> config::Seed;
326 break;
327
328 case 'R': // MaxOuterStep
329 cout << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ";
330 cin >> config::MaxOuterStep;
331 break;
332 case 'T': // OutVisStep
333 cout << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ";
334 cin >> config::OutVisStep;
335 break;
336 case 'U': // OutSrcStep
337 cout << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ";
338 cin >> config::OutSrcStep;
339 break;
340 case 'X': // MaxPsiStep
341 cout << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ";
342 cin >> config::MaxPsiStep;
343 break;
344 case 'Y': // EpsWannier
345 cout << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ";
346 cin >> config::EpsWannier;
347 break;
348
349 case 'Z': // MaxMinStep
350 cout << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ";
351 cin >> config::MaxMinStep;
352 break;
353 case 'a': // RelEpsTotalEnergy
354 cout << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ";
355 cin >> config::RelEpsTotalEnergy;
356 break;
357 case 'b': // RelEpsKineticEnergy
358 cout << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ";
359 cin >> config::RelEpsKineticEnergy;
360 break;
361 case 'c': // MaxMinStopStep
362 cout << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ";
363 cin >> config::MaxMinStopStep;
364 break;
365 case 'e': // MaxInitMinStep
366 cout << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ";
367 cin >> config::MaxInitMinStep;
368 break;
369 case 'f': // InitRelEpsTotalEnergy
370 cout << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ";
371 cin >> config::InitRelEpsTotalEnergy;
372 break;
373 case 'g': // InitRelEpsKineticEnergy
374 cout << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ";
375 cin >> config::InitRelEpsKineticEnergy;
376 break;
377 case 'h': // InitMaxMinStopStep
378 cout << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ";
379 cin >> config::InitMaxMinStopStep;
380 break;
381
382 case 'j': // BoxLength
383 cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
384 for (int i=0;i<6;i++) {
385 cout << Verbose(0) << "Cell size" << i << ": ";
386 cin >> mol->cell_size[i];
387 }
388 break;
389
390 case 'k': // ECut
391 cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
392 cin >> config::ECut;
393 break;
394 case 'l': // MaxLevel
395 cout << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ";
396 cin >> config::MaxLevel;
397 break;
398 case 'm': // RiemannTensor
399 cout << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ";
400 cin >> config::RiemannTensor;
401 break;
402 case 'n': // LevRFactor
403 cout << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ";
404 cin >> config::LevRFactor;
405 break;
406 case 'o': // RiemannLevel
407 cout << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ";
408 cin >> config::RiemannLevel;
409 break;
410 case 'p': // Lev0Factor
411 cout << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ";
412 cin >> config::Lev0Factor;
413 break;
414 case 'r': // RTActualUse
415 cout << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ";
416 cin >> config::RTActualUse;
417 break;
418 case 's': // PsiType
419 cout << Verbose(0) << "Old: " << config::PsiType << "\t new: ";
420 cin >> config::PsiType;
421 break;
422 case 't': // MaxPsiDouble
423 cout << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ";
424 cin >> config::MaxPsiDouble;
425 break;
426 case 'u': // PsiMaxNoUp
427 cout << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ";
428 cin >> config::PsiMaxNoUp;
429 break;
430 case 'v': // PsiMaxNoDown
431 cout << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ";
432 cin >> config::PsiMaxNoDown;
433 break;
434 case 'w': // AddPsis
435 cout << Verbose(0) << "Old: " << config::AddPsis << "\t new: ";
436 cin >> config::AddPsis;
437 break;
438
439 case 'x': // RCut
440 cout << Verbose(0) << "Old: " << config::RCut << "\t new: ";
441 cin >> config::RCut;
442 break;
443 case 'y': // StructOpt
444 cout << Verbose(0) << "Old: " << config::StructOpt << "\t new: ";
445 cin >> config::StructOpt;
446 break;
447 case 'z': // IsAngstroem
448 cout << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ";
449 cin >> config::IsAngstroem;
450 break;
451 case 'i': // RelativeCoord
452 cout << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ";
453 cin >> config::RelativeCoord;
454 break;
455 };
456 } while (choice != 'q');
457};
458
459/** Tests whether a given configuration file adhears to old or new syntax.
460 * \param *filename filename of config file to be tested
461 * \param *periode pointer to a periodentafel class with all elements
462 * \param *mol pointer to molecule containing all atoms of the molecule
463 * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax
464 */
465int config::TestSyntax(char *filename, periodentafel *periode, molecule *mol)
466{
467 int test;
468 ifstream file(filename);
469
470 // search file for keyword: ProcPEGamma (new syntax)
471 if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
472 file.close();
473 return 1;
474 }
475 // search file for keyword: ProcsGammaPsi (old syntax)
476 if (ParseForParameter(1,&file,"ProcsGammaPsi", 0, 1, 1, int_type, &test, 1, optional)) {
477 file.close();
478 return 0;
479 }
480 file.close();
481 return -1;
482}
483
484/** Returns private config::IsAngstroem.
485 * \return IsAngstroem
486 */
487bool config::GetIsAngstroem() const
488{
489 return (IsAngstroem == 1);
490};
491
492/** Returns private config::*defaultpath.
493 * \return *defaultpath
494 */
495char * config::GetDefaultPath() const
496{
497 return defaultpath;
498};
499
500
501/** Returns private config::*defaultpath.
502 * \return *defaultpath
503 */
504void config::SetDefaultPath(const char *path)
505{
506 strcpy(defaultpath, path);
507};
508
509/** Retrieves the path in the given config file name.
510 * \param filename config file string
511 */
512void config::RetrieveConfigPathAndName(string filename)
513{
514 char *ptr = NULL;
515 char *buffer = new char[MAXSTRINGSIZE];
516 strncpy(buffer, filename.c_str(), MAXSTRINGSIZE);
517 int last = -1;
518 for(last=MAXSTRINGSIZE;last--;) {
519 if (buffer[last] == '/')
520 break;
521 }
522 if (last == -1) { // no path in front, set to local directory.
523 strcpy(configpath, "./");
524 ptr = buffer;
525 } else {
526 strncpy(configpath, buffer, last+1);
527 ptr = &buffer[last+1];
528 if (last < 254)
529 configpath[last+1]='\0';
530 }
531 strcpy(configname, ptr);
532 cout << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl;
533 delete[](buffer);
534};
535
536
537/** Initializes config file structure by loading elements from a give file.
538 * \param *file input file stream being the opened config file
539 * \param *periode pointer to a periodentafel class with all elements
540 * \param *mol pointer to molecule containing all atoms of the molecule
541 */
542void config::Load(char *filename, periodentafel *periode, molecule *mol)
543{
544 ifstream *file = new ifstream(filename);
545 if (file == NULL) {
546 cerr << "ERROR: config file " << filename << " missing!" << endl;
547 return;
548 }
549 RetrieveConfigPathAndName(filename);
550 // ParseParameters
551
552 /* Oeffne Hauptparameterdatei */
553 int di;
554 double BoxLength[9];
555 string zeile;
556 string dummy;
557 element *elementhash[MAX_ELEMENTS];
558 char name[MAX_ELEMENTS];
559 char keyword[MAX_ELEMENTS];
560 int Z, No[MAX_ELEMENTS];
561 int verbose = 0;
562 double value[3];
563
564 InitThermostats(file);
565
566 /* Namen einlesen */
567
568 ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
569 ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
570 ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
571 ParseForParameter(verbose,file,"ProcPEGamma", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
572 ParseForParameter(verbose,file,"ProcPEPsi", 0, 1, 1, int_type, &(config::ProcPEPsi), 1, critical);
573
574 if (!ParseForParameter(verbose,file,"Seed", 0, 1, 1, int_type, &(config::Seed), 1, optional))
575 config::Seed = 1;
576
577 if(!ParseForParameter(verbose,file,"DoOutOrbitals", 0, 1, 1, int_type, &(config::DoOutOrbitals), 1, optional)) {
578 config::DoOutOrbitals = 0;
579 } else {
580 if (config::DoOutOrbitals < 0) config::DoOutOrbitals = 0;
581 if (config::DoOutOrbitals > 1) config::DoOutOrbitals = 1;
582 }
583 ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
584 if (config::DoOutVis < 0) config::DoOutVis = 0;
585 if (config::DoOutVis > 1) config::DoOutVis = 1;
586 if (!ParseForParameter(verbose,file,"VectorPlane", 0, 1, 1, int_type, &(config::VectorPlane), 1, optional))
587 config::VectorPlane = -1;
588 if (!ParseForParameter(verbose,file,"VectorCut", 0, 1, 1, double_type, &(config::VectorCut), 1, optional))
589 config::VectorCut = 0.;
590 ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
591 if (config::DoOutMes < 0) config::DoOutMes = 0;
592 if (config::DoOutMes > 1) config::DoOutMes = 1;
593 if (!ParseForParameter(verbose,file,"DoOutCurr", 0, 1, 1, int_type, &(config::DoOutCurrent), 1, optional))
594 config::DoOutCurrent = 0;
595 if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
596 if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
597 ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
598 if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
599 if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
600 if(!ParseForParameter(verbose,file,"DoWannier", 0, 1, 1, int_type, &(config::DoWannier), 1, optional)) {
601 config::DoWannier = 0;
602 } else {
603 if (config::DoWannier < 0) config::DoWannier = 0;
604 if (config::DoWannier > 1) config::DoWannier = 1;
605 }
606 if(!ParseForParameter(verbose,file,"CommonWannier", 0, 1, 1, int_type, &(config::CommonWannier), 1, optional)) {
607 config::CommonWannier = 0;
608 } else {
609 if (config::CommonWannier < 0) config::CommonWannier = 0;
610 if (config::CommonWannier > 4) config::CommonWannier = 4;
611 }
612 if(!ParseForParameter(verbose,file,"SawtoothStart", 0, 1, 1, double_type, &(config::SawtoothStart), 1, optional)) {
613 config::SawtoothStart = 0.01;
614 } else {
615 if (config::SawtoothStart < 0.) config::SawtoothStart = 0.;
616 if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
617 }
618
619 if (ParseForParameter(verbose,file,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))
620 if (config::DoConstrainedMD < 0)
621 config::DoConstrainedMD = 0;
622 ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
623 if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
624 config::Deltat = 1;
625 ParseForParameter(verbose,file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
626 ParseForParameter(verbose,file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
627 ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
628 //ParseForParameter(verbose,file,"Thermostat", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
629 if (!ParseForParameter(verbose,file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
630 config::EpsWannier = 1e-8;
631
632 // stop conditions
633 //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
634 ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
635 if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
636
637 ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
638 ParseForParameter(verbose,file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
639 ParseForParameter(verbose,file,"RelEpsKineticE", 0, 1, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
640 ParseForParameter(verbose,file,"MaxMinStopStep", 0, 1, 1, int_type, &(config::MaxMinStopStep), 1, critical);
641 ParseForParameter(verbose,file,"MaxMinGapStopStep", 0, 1, 1, int_type, &(config::MaxMinGapStopStep), 1, critical);
642 if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
643 if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
644 if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
645
646 ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
647 ParseForParameter(verbose,file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
648 ParseForParameter(verbose,file,"InitRelEpsKineticE", 0, 1, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
649 ParseForParameter(verbose,file,"InitMaxMinStopStep", 0, 1, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
650 ParseForParameter(verbose,file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &(config::InitMaxMinGapStopStep), 1, critical);
651 if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
652 if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
653 if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
654
655 // Unit cell and magnetic field
656 ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
657 mol->cell_size[0] = BoxLength[0];
658 mol->cell_size[1] = BoxLength[3];
659 mol->cell_size[2] = BoxLength[4];
660 mol->cell_size[3] = BoxLength[6];
661 mol->cell_size[4] = BoxLength[7];
662 mol->cell_size[5] = BoxLength[8];
663 if (1) fprintf(stderr,"\n");
664
665 ParseForParameter(verbose,file,"DoPerturbation", 0, 1, 1, int_type, &(config::DoPerturbation), 1, optional);
666 ParseForParameter(verbose,file,"DoOutNICS", 0, 1, 1, int_type, &(config::DoOutNICS), 1, optional);
667 if (!ParseForParameter(verbose,file,"DoFullCurrent", 0, 1, 1, int_type, &(config::DoFullCurrent), 1, optional))
668 config::DoFullCurrent = 0;
669 if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
670 if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
671 if (config::DoOutNICS < 0) config::DoOutNICS = 0;
672 if (config::DoOutNICS > 2) config::DoOutNICS = 2;
673 if (config::DoPerturbation == 0) {
674 config::DoFullCurrent = 0;
675 config::DoOutNICS = 0;
676 }
677
678 ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
679 ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
680 ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
681 if (config::Lev0Factor < 2) {
682 config::Lev0Factor = 2;
683 }
684 ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
685 if (di >= 0 && di < 2) {
686 config::RiemannTensor = di;
687 } else {
688 fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
689 exit(1);
690 }
691 switch (config::RiemannTensor) {
692 case 0: //UseNoRT
693 if (config::MaxLevel < 2) {
694 config::MaxLevel = 2;
695 }
696 config::LevRFactor = 2;
697 config::RTActualUse = 0;
698 break;
699 case 1: // UseRT
700 if (config::MaxLevel < 3) {
701 config::MaxLevel = 3;
702 }
703 ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
704 if (config::RiemannLevel < 2) {
705 config::RiemannLevel = 2;
706 }
707 if (config::RiemannLevel > config::MaxLevel-1) {
708 config::RiemannLevel = config::MaxLevel-1;
709 }
710 ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
711 if (config::LevRFactor < 2) {
712 config::LevRFactor = 2;
713 }
714 config::Lev0Factor = 2;
715 config::RTActualUse = 2;
716 break;
717 }
718 ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
719 if (di >= 0 && di < 2) {
720 config::PsiType = di;
721 } else {
722 fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
723 exit(1);
724 }
725 switch (config::PsiType) {
726 case 0: // SpinDouble
727 ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
728 ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
729 break;
730 case 1: // SpinUpDown
731 if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
732 ParseForParameter(verbose,file,"PsiMaxNoUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
733 ParseForParameter(verbose,file,"PsiMaxNoDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
734 ParseForParameter(verbose,file,"AddPsis", 0, 1, 1, int_type, &(config::AddPsis), 1, optional);
735 break;
736 }
737
738 // IonsInitRead
739
740 ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
741 ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
742 ParseForParameter(verbose,file,"MaxTypes", 0, 1, 1, int_type, &(config::MaxTypes), 1, critical);
743 if (!ParseForParameter(verbose,file,"RelativeCoord", 0, 1, 1, int_type, &(config::RelativeCoord) , 1, optional))
744 config::RelativeCoord = 0;
745 if (!ParseForParameter(verbose,file,"StructOpt", 0, 1, 1, int_type, &(config::StructOpt), 1, optional))
746 config::StructOpt = 0;
747 if (MaxTypes == 0) {
748 cerr << "There are no atoms according to MaxTypes in this config file." << endl;
749 } else {
750 // prescan number of ions per type
751 cout << Verbose(0) << "Prescanning ions per type: " << endl;
752 for (int i=0; i < config::MaxTypes; i++) {
753 sprintf(name,"Ion_Type%i",i+1);
754 ParseForParameter(verbose,file, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
755 ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
756 elementhash[i] = periode->FindElement(Z);
757 cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
758 }
759 int repetition = 0; // which repeated keyword shall be read
760
761 map<int, atom *> AtomList[config::MaxTypes];
762 if (!FastParsing) {
763 // parse in trajectories
764 bool status = true;
765 atom *neues = NULL;
766 while (status) {
767 cout << "Currently parsing MD step " << repetition << "." << endl;
768 for (int i=0; i < config::MaxTypes; i++) {
769 sprintf(name,"Ion_Type%i",i+1);
770 for(int j=0;j<No[i];j++) {
771 sprintf(keyword,"%s_%i",name, j+1);
772 if (repetition == 0) {
773 neues = new atom();
774 AtomList[i][j] = neues;
775 neues->type = elementhash[i]; // find element type
776 mol->AddAtom(neues);
777 } else
778 neues = AtomList[i][j];
779 status = (status &&
780 ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
781 ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
782 ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
783 ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
784 if (!status) break;
785
786 // check size of vectors
787 if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
788 //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
789 mol->Trajectories[neues].R.resize(repetition+10);
790 mol->Trajectories[neues].U.resize(repetition+10);
791 mol->Trajectories[neues].F.resize(repetition+10);
792 }
793
794 // put into trajectories list
795 for (int d=0;d<NDIM;d++)
796 mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
797
798 // parse velocities if present
799 if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
800 neues->v.x[0] = 0.;
801 if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
802 neues->v.x[1] = 0.;
803 if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
804 neues->v.x[2] = 0.;
805 for (int d=0;d<NDIM;d++)
806 mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
807
808 // parse forces if present
809 if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
810 value[0] = 0.;
811 if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
812 value[1] = 0.;
813 if(!ParseForParameter(verbose,file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
814 value[2] = 0.;
815 for (int d=0;d<NDIM;d++)
816 mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
817
818 // cout << "Parsed position of step " << (repetition) << ": (";
819 // for (int d=0;d<NDIM;d++)
820 // cout << mol->Trajectories[neues].R.at(repetition).x[d] << " "; // next step
821 // cout << ")\t(";
822 // for (int d=0;d<NDIM;d++)
823 // cout << mol->Trajectories[neues].U.at(repetition).x[d] << " "; // next step
824 // cout << ")\t(";
825 // for (int d=0;d<NDIM;d++)
826 // cout << mol->Trajectories[neues].F.at(repetition).x[d] << " "; // next step
827 // cout << ")" << endl;
828 }
829 }
830 repetition++;
831 }
832 repetition--;
833 cout << "Found " << repetition << " trajectory steps." << endl;
834 mol->MDSteps = repetition;
835 } else {
836 // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
837 repetition = 0;
838 while ( ParseForParameter(verbose,file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
839 ParseForParameter(verbose,file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
840 ParseForParameter(verbose,file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
841 repetition++;
842 cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
843 // parse in molecule coordinates
844 for (int i=0; i < config::MaxTypes; i++) {
845 sprintf(name,"Ion_Type%i",i+1);
846 for(int j=0;j<No[i];j++) {
847 sprintf(keyword,"%s_%i",name, j+1);
848 atom *neues = new atom();
849 // then parse for each atom the coordinates as often as present
850 ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
851 ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
852 ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
853 ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
854 if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
855 neues->v.x[0] = 0.;
856 if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
857 neues->v.x[1] = 0.;
858 if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
859 neues->v.x[2] = 0.;
860 // here we don't care if forces are present (last in trajectories is always equal to current position)
861 neues->type = elementhash[i]; // find element type
862 mol->AddAtom(neues);
863 }
864 }
865 }
866 }
867 file->close();
868 delete(file);
869};
870
871/** Initializes config file structure by loading elements from a give file with old pcp syntax.
872 * \param *file input file stream being the opened config file with old pcp syntax
873 * \param *periode pointer to a periodentafel class with all elements
874 * \param *mol pointer to molecule containing all atoms of the molecule
875 */
876void config::LoadOld(char *filename, periodentafel *periode, molecule *mol)
877{
878 ifstream *file = new ifstream(filename);
879 if (file == NULL) {
880 cerr << "ERROR: config file " << filename << " missing!" << endl;
881 return;
882 }
883 RetrieveConfigPathAndName(filename);
884 // ParseParameters
885
886 /* Oeffne Hauptparameterdatei */
887 int l, i, di;
888 double a,b;
889 double BoxLength[9];
890 string zeile;
891 string dummy;
892 element *elementhash[128];
893 int Z, No, AtomNo, found;
894 int verbose = 0;
895
896 /* Namen einlesen */
897
898 ParseForParameter(verbose,file, "mainname", 0, 1, 1, string_type, (config::mainname), 1, critical);
899 ParseForParameter(verbose,file, "defaultpath", 0, 1, 1, string_type, (config::defaultpath), 1, critical);
900 ParseForParameter(verbose,file, "pseudopotpath", 0, 1, 1, string_type, (config::pseudopotpath), 1, critical);
901 ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 1, 1, int_type, &(config::ProcPEGamma), 1, critical);
902 ParseForParameter(verbose,file, "ProcsGammaPsi", 0, 2, 1, int_type, &(config::ProcPEPsi), 1, critical);
903 config::Seed = 1;
904 config::DoOutOrbitals = 0;
905 ParseForParameter(verbose,file,"DoOutVis", 0, 1, 1, int_type, &(config::DoOutVis), 1, critical);
906 if (config::DoOutVis < 0) config::DoOutVis = 0;
907 if (config::DoOutVis > 1) config::DoOutVis = 1;
908 config::VectorPlane = -1;
909 config::VectorCut = 0.;
910 ParseForParameter(verbose,file,"DoOutMes", 0, 1, 1, int_type, &(config::DoOutMes), 1, critical);
911 if (config::DoOutMes < 0) config::DoOutMes = 0;
912 if (config::DoOutMes > 1) config::DoOutMes = 1;
913 config::DoOutCurrent = 0;
914 ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
915 if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
916 if (config::UseAddGramSch > 2) config::UseAddGramSch = 2;
917 config::CommonWannier = 0;
918 config::SawtoothStart = 0.01;
919
920 ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, double_type, &(config::MaxOuterStep), 1, critical);
921 ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
922 ParseForParameter(verbose,file,"VisOuterStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
923 ParseForParameter(verbose,file,"VisSrcOuterStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
924 ParseForParameter(verbose,file,"TargetTemp", 0, 1, 1, double_type, &(config::TargetTemp), 1, optional);
925 ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
926 config::EpsWannier = 1e-8;
927
928 // stop conditions
929 //if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
930 ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
931 if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
932
933 ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
934 ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
935 ParseForParameter(verbose,file,"MaxMinStep", 0, 3, 1, double_type, &(config::RelEpsKineticEnergy), 1, critical);
936 ParseForParameter(verbose,file,"MaxMinStep", 0, 4, 1, int_type, &(config::MaxMinStopStep), 1, critical);
937 if (config::MaxMinStep <= 0) config::MaxMinStep = config::MaxPsiStep;
938 if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
939 config::MaxMinGapStopStep = 1;
940
941 ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
942 ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
943 ParseForParameter(verbose,file,"MaxInitMinStep", 0, 3, 1, double_type, &(config::InitRelEpsKineticEnergy), 1, critical);
944 ParseForParameter(verbose,file,"MaxInitMinStep", 0, 4, 1, int_type, &(config::InitMaxMinStopStep), 1, critical);
945 if (config::MaxInitMinStep <= 0) config::MaxInitMinStep = config::MaxPsiStep;
946 if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
947 config::InitMaxMinGapStopStep = 1;
948
949 ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
950 mol->cell_size[0] = BoxLength[0];
951 mol->cell_size[1] = BoxLength[3];
952 mol->cell_size[2] = BoxLength[4];
953 mol->cell_size[3] = BoxLength[6];
954 mol->cell_size[4] = BoxLength[7];
955 mol->cell_size[5] = BoxLength[8];
956 if (1) fprintf(stderr,"\n");
957 config::DoPerturbation = 0;
958 config::DoFullCurrent = 0;
959
960 ParseForParameter(verbose,file,"ECut", 0, 1, 1, double_type, &(config::ECut), 1, critical);
961 ParseForParameter(verbose,file,"MaxLevel", 0, 1, 1, int_type, &(config::MaxLevel), 1, critical);
962 ParseForParameter(verbose,file,"Level0Factor", 0, 1, 1, int_type, &(config::Lev0Factor), 1, critical);
963 if (config::Lev0Factor < 2) {
964 config::Lev0Factor = 2;
965 }
966 ParseForParameter(verbose,file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
967 if (di >= 0 && di < 2) {
968 config::RiemannTensor = di;
969 } else {
970 fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
971 exit(1);
972 }
973 switch (config::RiemannTensor) {
974 case 0: //UseNoRT
975 if (config::MaxLevel < 2) {
976 config::MaxLevel = 2;
977 }
978 config::LevRFactor = 2;
979 config::RTActualUse = 0;
980 break;
981 case 1: // UseRT
982 if (config::MaxLevel < 3) {
983 config::MaxLevel = 3;
984 }
985 ParseForParameter(verbose,file,"RiemannLevel", 0, 1, 1, int_type, &(config::RiemannLevel), 1, critical);
986 if (config::RiemannLevel < 2) {
987 config::RiemannLevel = 2;
988 }
989 if (config::RiemannLevel > config::MaxLevel-1) {
990 config::RiemannLevel = config::MaxLevel-1;
991 }
992 ParseForParameter(verbose,file,"LevRFactor", 0, 1, 1, int_type, &(config::LevRFactor), 1, critical);
993 if (config::LevRFactor < 2) {
994 config::LevRFactor = 2;
995 }
996 config::Lev0Factor = 2;
997 config::RTActualUse = 2;
998 break;
999 }
1000 ParseForParameter(verbose,file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
1001 if (di >= 0 && di < 2) {
1002 config::PsiType = di;
1003 } else {
1004 fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
1005 exit(1);
1006 }
1007 switch (config::PsiType) {
1008 case 0: // SpinDouble
1009 ParseForParameter(verbose,file,"MaxPsiDouble", 0, 1, 1, int_type, &(config::MaxPsiDouble), 1, critical);
1010 config::AddPsis = 0;
1011 break;
1012 case 1: // SpinUpDown
1013 if (config::ProcPEGamma % 2) config::ProcPEGamma*=2;
1014 ParseForParameter(verbose,file,"MaxPsiUp", 0, 1, 1, int_type, &(config::PsiMaxNoUp), 1, critical);
1015 ParseForParameter(verbose,file,"MaxPsiDown", 0, 1, 1, int_type, &(config::PsiMaxNoDown), 1, critical);
1016 config::AddPsis = 0;
1017 break;
1018 }
1019
1020 // IonsInitRead
1021
1022 ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
1023 ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
1024 config::RelativeCoord = 0;
1025 config::StructOpt = 0;
1026
1027 // Routine from builder.cpp
1028
1029
1030 for (i=MAX_ELEMENTS;i--;)
1031 elementhash[i] = NULL;
1032 cout << Verbose(0) << "Parsing Ions ..." << endl;
1033 No=0;
1034 found = 0;
1035 while (getline(*file,zeile,'\n')) {
1036 if (zeile.find("Ions_Data") == 0) {
1037 cout << Verbose(1) << "found Ions_Data...begin parsing" << endl;
1038 found ++;
1039 }
1040 if (found > 0) {
1041 if (zeile.find("Ions_Data") == 0)
1042 getline(*file,zeile,'\n'); // read next line and parse this one
1043 istringstream input(zeile);
1044 input >> AtomNo; // number of atoms
1045 input >> Z; // atomic number
1046 input >> a;
1047 input >> l;
1048 input >> l;
1049 input >> b; // element mass
1050 elementhash[No] = periode->FindElement(Z);
1051 cout << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:" << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl;
1052 for(i=0;i<AtomNo;i++) {
1053 if (!getline(*file,zeile,'\n')) {// parse on and on
1054 cout << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl;
1055 // return 1;
1056 } else {
1057 //cout << Verbose(2) << "Reading line: " << zeile << endl;
1058 }
1059 istringstream input2(zeile);
1060 atom *neues = new atom();
1061 input2 >> neues->x.x[0]; // x
1062 input2 >> neues->x.x[1]; // y
1063 input2 >> neues->x.x[2]; // z
1064 input2 >> l;
1065 neues->type = elementhash[No]; // find element type
1066 mol->AddAtom(neues);
1067 }
1068 No++;
1069 }
1070 }
1071 file->close();
1072 delete(file);
1073};
1074
1075/** Stores all elements of config structure from which they can be re-read.
1076 * \param output open output *file stream to write to
1077 * \param *periode pointer to a periodentafel class with all elements
1078 * \param *mol pointer to molecule containing all atoms of the molecule
1079 */
1080bool config::Save(ofstream *output, periodentafel *periode, molecule *mol) const
1081{
1082 bool result = true;
1083 // bring MaxTypes up to date
1084 mol->CountElements();
1085 if (output != NULL) {
1086 *output << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
1087 *output << endl;
1088 *output << "mainname\t" << config::mainname << "\t# programm name (for runtime files)" << endl;
1089 *output << "defaultpath\t" << config::defaultpath << "\t# where to put files during runtime" << endl;
1090 *output << "pseudopotpath\t" << config::pseudopotpath << "\t# where to find pseudopotentials" << endl;
1091 *output << endl;
1092 *output << "ProcPEGamma\t" << config::ProcPEGamma << "\t# for parallel computing: share constants" << endl;
1093 *output << "ProcPEPsi\t" << config::ProcPEPsi << "\t# for parallel computing: share wave functions" << endl;
1094 *output << "DoOutVis\t" << config::DoOutVis << "\t# Output data for OpenDX" << endl;
1095 *output << "DoOutMes\t" << config::DoOutMes << "\t# Output data for measurements" << endl;
1096 *output << "DoOutOrbitals\t" << config::DoOutOrbitals << "\t# Output all Orbitals" << endl;
1097 *output << "DoOutCurr\t" << config::DoOutCurrent << "\t# Ouput current density for OpenDx" << endl;
1098 *output << "DoOutNICS\t" << config::DoOutNICS << "\t# Output Nucleus independent current shieldings" << endl;
1099 *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
1100 *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
1101 *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
1102 *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
1103 switch(Thermostat) {
1104 default:
1105 case None:
1106 break;
1107 case Woodcock:
1108 *output << ScaleTempStep;
1109 break;
1110 case Gaussian:
1111 *output << ScaleTempStep;
1112 break;
1113 case Langevin:
1114 *output << TempFrequency << "\t" << alpha;
1115 break;
1116 case Berendsen:
1117 *output << TempFrequency;
1118 break;
1119 case NoseHoover:
1120 *output << HooverMass;
1121 break;
1122 };
1123 *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
1124 *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
1125 *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
1126 *output << "VectorPlane\t" << config::VectorPlane << "\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot" << endl;
1127 *output << "VectorCut\t" << config::VectorCut << "\t# Cut plane axis value" << endl;
1128 *output << "AddGramSch\t" << config::UseAddGramSch << "\t# Additional GramSchmidtOrtogonalization to be safe" << endl;
1129 *output << "Seed\t\t" << config::Seed << "\t# initial value for random seed for Psi coefficients" << endl;
1130 *output << endl;
1131 *output << "MaxOuterStep\t" << config::MaxOuterStep << "\t# number of MolecularDynamics/Structure optimization steps" << endl;
1132 *output << "Deltat\t" << config::Deltat << "\t# time per MD step" << endl;
1133 *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
1134 *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
1135 *output << "TargetTemp\t" << config::TargetTemp << "\t# Target temperature" << endl;
1136 *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
1137 *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
1138 *output << endl;
1139 *output << "# Values specifying when to stop" << endl;
1140 *output << "MaxMinStep\t" << config::MaxMinStep << "\t# Maximum number of steps" << endl;
1141 *output << "RelEpsTotalE\t" << config::RelEpsTotalEnergy << "\t# relative change in total energy" << endl;
1142 *output << "RelEpsKineticE\t" << config::RelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
1143 *output << "MaxMinStopStep\t" << config::MaxMinStopStep << "\t# check every ..th steps" << endl;
1144 *output << "MaxMinGapStopStep\t" << config::MaxMinGapStopStep << "\t# check every ..th steps" << endl;
1145 *output << endl;
1146 *output << "# Values specifying when to stop for INIT, otherwise same as above" << endl;
1147 *output << "MaxInitMinStep\t" << config::MaxInitMinStep << "\t# Maximum number of steps" << endl;
1148 *output << "InitRelEpsTotalE\t" << config::InitRelEpsTotalEnergy << "\t# relative change in total energy" << endl;
1149 *output << "InitRelEpsKineticE\t" << config::InitRelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
1150 *output << "InitMaxMinStopStep\t" << config::InitMaxMinStopStep << "\t# check every ..th steps" << endl;
1151 *output << "InitMaxMinGapStopStep\t" << config::InitMaxMinGapStopStep << "\t# check every ..th steps" << endl;
1152 *output << endl;
1153 *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
1154 *output << mol->cell_size[0] << "\t" << endl;
1155 *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
1156 *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
1157 // FIXME
1158 *output << endl;
1159 *output << "ECut\t\t" << config::ECut << "\t# energy cutoff for discretization in Hartrees" << endl;
1160 *output << "MaxLevel\t" << config::MaxLevel << "\t# number of different levels in the code, >=2" << endl;
1161 *output << "Level0Factor\t" << config::Lev0Factor << "\t# factor by which node number increases from S to 0 level" << endl;
1162 *output << "RiemannTensor\t" << config::RiemannTensor << "\t# (Use metric)" << endl;
1163 switch (config::RiemannTensor) {
1164 case 0: //UseNoRT
1165 break;
1166 case 1: // UseRT
1167 *output << "RiemannLevel\t" << config::RiemannLevel << "\t# Number of Riemann Levels" << endl;
1168 *output << "LevRFactor\t" << config::LevRFactor << "\t# factor by which node number increases from 0 to R level from" << endl;
1169 break;
1170 }
1171 *output << "PsiType\t\t" << config::PsiType << "\t# 0 - doubly occupied, 1 - SpinUp,SpinDown" << endl;
1172 // write out both types for easier changing afterwards
1173 // switch (PsiType) {
1174 // case 0:
1175 *output << "MaxPsiDouble\t" << config::MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
1176 // break;
1177 // case 1:
1178 *output << "PsiMaxNoUp\t" << config::PsiMaxNoUp << "\t# here: specifying maximum number of SpinUp-states" << endl;
1179 *output << "PsiMaxNoDown\t" << config::PsiMaxNoDown << "\t# here: specifying maximum number of SpinDown-states" << endl;
1180 // break;
1181 // }
1182 *output << "AddPsis\t\t" << config::AddPsis << "\t# Additional unoccupied Psis for bandgap determination" << endl;
1183 *output << endl;
1184 *output << "RCut\t\t" << config::RCut << "\t# R-cut for the ewald summation" << endl;
1185 *output << "StructOpt\t" << config::StructOpt << "\t# Do structure optimization beforehand" << endl;
1186 *output << "IsAngstroem\t" << config::IsAngstroem << "\t# 0 - Bohr, 1 - Angstroem" << endl;
1187 *output << "RelativeCoord\t" << config::RelativeCoord << "\t# whether ion coordinates are relative (1) or absolute (0)" << endl;
1188 *output << "MaxTypes\t" << mol->ElementCount << "\t# maximum number of different ion types" << endl;
1189 *output << endl;
1190 result = result && mol->Checkout(output);
1191 if (mol->MDSteps <=1 )
1192 result = result && mol->Output(output);
1193 else
1194 result = result && mol->OutputTrajectories(output);
1195 return result;
1196 } else
1197 return false;
1198};
1199
1200/** Stores all elements in a MPQC input file.
1201 * Note that this format cannot be parsed again.
1202 * \param output open output *file stream to write to
1203 * \param *mol pointer to molecule containing all atoms of the molecule
1204 */
1205bool config::SaveMPQC(ofstream *output, molecule *mol) const
1206{
1207 int ElementNo = 0;
1208 int AtomNo;
1209 atom *Walker = NULL;
1210 element *runner = mol->elemente->start;
1211 Vector *center = NULL;
1212
1213 *output << "% Created by MoleCuilder" << endl;
1214 *output << "mpqc: (" << endl;
1215 *output << "\tsavestate = no" << endl;
1216 *output << "\tdo_gradient = yes" << endl;
1217 *output << "\tmole<CLHF>: (" << endl;
1218 *output << "\t\tmolecule<Molecule>: (" << endl;
1219 *output << "\t\t\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
1220 *output << "\t\t\t{ atoms geometry } = {" << endl;
1221 center = mol->DetermineCenterOfAll(output);
1222 // output of atoms
1223 while (runner->next != mol->elemente->end) { // go through every element
1224 runner = runner->next;
1225 if (mol->ElementsInMolecule[runner->Z]) { // if this element got atoms
1226 ElementNo++;
1227 AtomNo = 0;
1228 Walker = mol->start;
1229 while (Walker->next != mol->end) { // go through every atom of this element
1230 Walker = Walker->next;
1231 if (Walker->type == runner) { // if this atom fits to element
1232 AtomNo++;
1233 *output << "\t\t\t\t" << Walker->type->symbol << " [ " << Walker->x.x[0]-center->x[0] << "\t" << Walker->x.x[1]-center->x[1] << "\t" << Walker->x.x[2]-center->x[2] << " ]" << endl;
1234 }
1235 }
1236 }
1237 }
1238 delete(center);
1239 *output << "\t\t\t}" << endl;
1240 *output << "\t\t)" << endl;
1241 *output << "\t\tbasis<GaussianBasisSet>: (" << endl;
1242 *output << "\t\t\tname = \"STO-3G\"" << endl;
1243 *output << "\t\t\tmolecule = $:mpqc:mole:molecule" << endl;
1244 *output << "\t\t)" << endl;
1245 *output << "\t)" << endl;
1246 *output << ")" << endl;
1247 return true;
1248};
1249
1250/** Reads parameter from a parsed file.
1251 * The file is either parsed for a certain keyword or if null is given for
1252 * the value in row yth and column xth. If the keyword was necessity#critical,
1253 * then an error is thrown and the programme aborted.
1254 * \warning value is modified (both in contents and position)!
1255 * \param verbose 1 - print found value to stderr, 0 - don't
1256 * \param file file to be parsed
1257 * \param name Name of value in file (at least 3 chars!)
1258 * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
1259 * (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
1260 * best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
1261 * counted from this unresetted position!)
1262 * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
1263 * \param yth In grid case specifying column number, otherwise the yth \a name matching line
1264 * \param type Type of the Parameter to be read
1265 * \param value address of the value to be read (must have been allocated)
1266 * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
1267 * \param critical necessity of this keyword being specified (optional, critical)
1268 * \return 1 - found, 0 - not found
1269 * \note Routine is taken from the pcp project and hack-a-slack adapted to C++
1270 */
1271int config::ParseForParameter(int verbose, ifstream *file, const char *name, int sequential, int const xth, int const yth, int type, void *value, int repetition, int critical) {
1272 int i,j; // loop variables
1273 int length = 0, maxlength = -1;
1274 long file_position = file->tellg(); // mark current position
1275 char *dummy1, *dummy, *free_dummy; // pointers in the line that is read in per step
1276 dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char), "config::ParseForParameter: *free_dummy");
1277
1278 //fprintf(stderr,"Parsing for %s\n",name);
1279 if (repetition == 0)
1280 //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
1281 return 0;
1282
1283 int line = 0; // marks line where parameter was found
1284 int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found
1285 while((found != repetition)) {
1286 dummy1 = dummy = free_dummy;
1287 do {
1288 file->getline(dummy1, 256); // Read the whole line
1289 if (file->eof()) {
1290 if ((critical) && (found == 0)) {
1291 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1292 //Error(InitReading, name);
1293 fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
1294 exit(255);
1295 } else {
1296 //if (!sequential)
1297 file->clear();
1298 file->seekg(file_position, ios::beg); // rewind to start position
1299 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1300 return 0;
1301 }
1302 }
1303 line++;
1304 } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
1305
1306 // C++ getline removes newline at end, thus re-add
1307 if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
1308 i = strlen(dummy1);
1309 dummy1[i] = '\n';
1310 dummy1[i+1] = '\0';
1311 }
1312 //fprintf(stderr,"line %i ends at %i, newline at %i\n", line, strlen(dummy1), strchr(dummy1,'\n')-free_dummy);
1313
1314 if (dummy1 == NULL) {
1315 if (verbose) fprintf(stderr,"Error reading line %i\n",line);
1316 } else {
1317 //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
1318 }
1319 // Seek for possible end of keyword on line if given ...
1320 if (name != NULL) {
1321 dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever's nearer
1322 if (dummy == NULL) {
1323 dummy = strchr(dummy1, ' '); // if not found seek for space
1324 while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary
1325 dummy++;
1326 }
1327 if (dummy == NULL) {
1328 dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
1329 //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
1330 //Free((void **)&free_dummy);
1331 //Error(FileOpenParams, NULL);
1332 } else {
1333 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
1334 }
1335 } else dummy = dummy1;
1336 // ... and check if it is the keyword!
1337 //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
1338 if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
1339 found++; // found the parameter!
1340 //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);
1341
1342 if (found == repetition) {
1343 for (i=0;i<xth;i++) { // i = rows
1344 if (type >= grid) {
1345 // grid structure means that grid starts on the next line, not right after keyword
1346 dummy1 = dummy = free_dummy;
1347 do {
1348 file->getline(dummy1, 256); // Read the whole line, skip commentary and empty ones
1349 if (file->eof()) {
1350 if ((critical) && (found == 0)) {
1351 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1352 //Error(InitReading, name);
1353 fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
1354 exit(255);
1355 } else {
1356 //if (!sequential)
1357 file->clear();
1358 file->seekg(file_position, ios::beg); // rewind to start position
1359 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1360 return 0;
1361 }
1362 }
1363 line++;
1364 } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
1365 if (dummy1 == NULL){
1366 if (verbose) fprintf(stderr,"Error reading line %i\n", line);
1367 } else {
1368 //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
1369 }
1370 } else { // simple int, strings or doubles start in the same line
1371 while ((*dummy == '\t') || (*dummy == ' ')) // skip interjacent tabs and spaces
1372 dummy++;
1373 }
1374 // C++ getline removes newline at end, thus re-add
1375 if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
1376 j = strlen(dummy1);
1377 dummy1[j] = '\n';
1378 dummy1[j+1] = '\0';
1379 }
1380
1381 int start = (type >= grid) ? 0 : yth-1 ;
1382 for (j=start;j<yth;j++) { // j = columns
1383 // check for lower triangular area and upper triangular area
1384 if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
1385 *((double *)value) = 0.0;
1386 fprintf(stderr,"%f\t",*((double *)value));
1387 value = (void *)((long)value + sizeof(double));
1388 //value += sizeof(double);
1389 } else {
1390 // otherwise we must skip all interjacent tabs and spaces and find next value
1391 dummy1 = dummy;
1392 dummy = strchr(dummy1, '\t'); // seek for tab or space
1393 if (dummy == NULL)
1394 dummy = strchr(dummy1, ' '); // if not found seek for space
1395 if (dummy == NULL) { // if still zero returned ...
1396 dummy = strchr(dummy1, '\n'); // ... at line end then
1397 if ((j < yth-1) && (type < 4)) { // check if xth value or not yet
1398 if (critical) {
1399 if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
1400 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1401 //return 0;
1402 exit(255);
1403 //Error(FileOpenParams, NULL);
1404 } else {
1405 //if (!sequential)
1406 file->clear();
1407 file->seekg(file_position, ios::beg); // rewind to start position
1408 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1409 return 0;
1410 }
1411 }
1412 } else {
1413 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
1414 }
1415 if (*dummy1 == '#') {
1416 // found comment, skipping rest of line
1417 //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
1418 if (!sequential) { // here we need it!
1419 file->seekg(file_position, ios::beg); // rewind to start position
1420 }
1421 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1422 return 0;
1423 }
1424 //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
1425 switch(type) {
1426 case (row_int):
1427 *((int *)value) = atoi(dummy1);
1428 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
1429 if (verbose) fprintf(stderr,"%i\t",*((int *)value));
1430 value = (void *)((long)value + sizeof(int));
1431 //value += sizeof(int);
1432 break;
1433 case(row_double):
1434 case(grid):
1435 case(lower_trigrid):
1436 case(upper_trigrid):
1437 *((double *)value) = atof(dummy1);
1438 if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
1439 if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
1440 value = (void *)((long)value + sizeof(double));
1441 //value += sizeof(double);
1442 break;
1443 case(double_type):
1444 *((double *)value) = atof(dummy1);
1445 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
1446 //value += sizeof(double);
1447 break;
1448 case(int_type):
1449 *((int *)value) = atoi(dummy1);
1450 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
1451 //value += sizeof(int);
1452 break;
1453 default:
1454 case(string_type):
1455 if (value != NULL) {
1456 //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
1457 maxlength = MAXSTRINGSIZE;
1458 length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
1459 strncpy((char *)value, dummy1, length); // copy as much
1460 ((char *)value)[length] = '\0'; // and set end marker
1461 if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
1462 //value += sizeof(char);
1463 } else {
1464 }
1465 break;
1466 }
1467 }
1468 while (*dummy == '\t')
1469 dummy++;
1470 }
1471 }
1472 }
1473 }
1474 }
1475 if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
1476 Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
1477 if (!sequential) {
1478 file->clear();
1479 file->seekg(file_position, ios::beg); // rewind to start position
1480 }
1481 //fprintf(stderr, "End of Parsing\n\n");
1482
1483 return (found); // true if found, false if not
1484}
Note: See TracBrowser for help on using the repository browser.