source: src/config.cpp@ d067d45

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 d067d45 was d067d45, checked in by Frederik Heber <heber@…>, 16 years ago

Merge branch 'MultipleMolecules'

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/atom.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/helpers.cpp
molecuilder/src/joiner.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp
molecuilder/src/vector.cpp
molecuilder/src/verbose.cpp

merges:

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