- Timestamp:
- Oct 19, 2008, 1:57:23 PM (16 years ago)
- Branches:
- 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
- Children:
- 84fc91
- Parents:
- f731ae
- Location:
- src
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified src/analyzer.cpp ¶
rf731ae reeec8f 25 25 periodentafel *periode = NULL; // and a period table of all elements 26 26 EnergyMatrix Energy; 27 EnergyMatrix EnergyFragments; 28 ForceMatrix Force; 29 ForceMatrix ForceFragments; 30 HessianMatrix Hessian; 31 HessianMatrix HessianFragments; 27 32 EnergyMatrix Hcorrection; 28 ForceMatrix Force;33 EnergyMatrix HcorrectionFragments; 29 34 ForceMatrix Shielding; 30 35 ForceMatrix ShieldingPAS; 31 36 EnergyMatrix Time; 32 EnergyMatrix EnergyFragments;33 EnergyMatrix HcorrectionFragments;34 ForceMatrix ForceFragments;35 37 ForceMatrix ShieldingFragments; 36 38 ForceMatrix ShieldingPASFragments; … … 49 51 stringstream yrange; 50 52 char *dir = NULL; 51 bool Hcorrected = true; 53 bool NoHCorrection = false; 54 bool NoHessian = false; 55 bool NoTime = false; 52 56 double norm; 53 57 int counter; … … 84 88 // ------------- Parse through all Fragment subdirs -------- 85 89 if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix,0,0)) return 1; 86 Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0); 90 if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX,0,0)) { 91 NoHCorrection = true; 92 cout << "No HCorrection file found, skipping these." << endl; 93 } 94 87 95 if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix,0,0)) return 1; 88 if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) return 1; 96 if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix,0,0)) { 97 NoHessian = true; 98 cout << "No Hessian file found, skipping these." << endl; 99 } 100 if (!Time.ParseFragmentMatrix(argv[1], dir, TimeSuffix, 10,1)) { 101 NoTime = true; 102 cout << "No speed file found, skipping these." << endl; 103 } 89 104 if (periode != NULL) { // also look for PAS values 90 105 if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1; … … 94 109 // ---------- Parse the TE Factors into an array ----------------- 95 110 if (!Energy.InitialiseIndices()) return 1; 96 if (Hcorrected) Hcorrection.InitialiseIndices(); 111 if (!NoHCorrection) 112 Hcorrection.InitialiseIndices(); 97 113 98 114 // ---------- Parse the Force indices into an array --------------- 99 115 if (!Force.ParseIndices(argv[1])) return 1; 100 116 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1; 101 if (!ForceFragments.ParseIndices(argv[1])) return 1; 117 if (!ForceFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1; 118 119 // ---------- Parse hessian indices into an array ----------------- 120 if (!NoHessian) { 121 if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1; 122 if (!HessianFragments.AllocateMatrix(Hessian.Header, Hessian.MatrixCounter, Hessian.RowCounter, Hessian.ColumnCounter)) return 1; 123 if (!HessianFragments.InitialiseIndices((class MatrixContainer *)&Force)) return 1; 124 } 102 125 103 126 // ---------- Parse the shielding indices into an array --------------- … … 117 140 // ---------- Parse fragment files created by 'joiner' into an array ------------- 118 141 if (!EnergyFragments.ParseFragmentMatrix(argv[1], dir, EnergyFragmentSuffix,0,0)) return 1; 119 if (Hcorrected) HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0); 142 if (!NoHCorrection) 143 HcorrectionFragments.ParseFragmentMatrix(argv[1], dir, HcorrectionFragmentSuffix,0,0); 120 144 if (!ForceFragments.ParseFragmentMatrix(argv[1], dir, ForceFragmentSuffix,0,0)) return 1; 145 if (!NoHessian) 146 if (!HessianFragments.ParseFragmentMatrix(argv[1], dir, HessianFragmentSuffix,0,0)) return 1; 121 147 if (periode != NULL) { // also look for PAS values 122 148 if (!ShieldingFragments.ParseFragmentMatrix(argv[1], dir, ShieldingFragmentSuffix, 1, 0)) return 1; … … 130 156 filename << argv[3] << "/" << "energy-forces.all"; 131 157 output.open(filename.str().c_str(), ios::out); 132 output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header << endl;158 output << endl << "Total Energy" << endl << "==============" << endl << Energy.Header[Energy.MatrixCounter] << endl; 133 159 for(int j=0;j<Energy.RowCounter[Energy.MatrixCounter];j++) { 134 160 for(int k=0;k<Energy.ColumnCounter[Energy.MatrixCounter];k++) … … 138 164 output << endl; 139 165 140 output << endl << "Total Forces" << endl << "===============" << endl << Force.Header << endl;166 output << endl << "Total Forces" << endl << "===============" << endl << Force.Header[Force.MatrixCounter] << endl; 141 167 for(int j=0;j<Force.RowCounter[Force.MatrixCounter];j++) { 142 168 for(int k=0;k<Force.ColumnCounter[Force.MatrixCounter];k++) … … 146 172 output << endl; 147 173 148 if (periode != NULL) { // also look for PAS values 149 output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header << endl; 174 if (!NoHessian) { 175 output << endl << "Total Hessian" << endl << "===============" << endl << Hessian.Header[Hessian.MatrixCounter] << endl; 176 for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) { 177 for(int k=0;k<Hessian.ColumnCounter[Hessian.MatrixCounter];k++) 178 output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t"; 179 output << endl; 180 } 181 output << endl; 182 } 183 184 if (periode != NULL) { // also look for PAS values 185 output << endl << "Total Shieldings" << endl << "===============" << endl << Shielding.Header[Hessian.MatrixCounter] << endl; 150 186 for(int j=0;j<Shielding.RowCounter[Shielding.MatrixCounter];j++) { 151 187 for(int k=0;k<Shielding.ColumnCounter[Shielding.MatrixCounter];k++) … … 155 191 output << endl; 156 192 157 output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header << endl;193 output << endl << "Total Shieldings PAS" << endl << "===============" << endl << ShieldingPAS.Header[ShieldingPAS.MatrixCounter] << endl; 158 194 for(int j=0;j<ShieldingPAS.RowCounter[ShieldingPAS.MatrixCounter];j++) { 159 195 for(int k=0;k<ShieldingPAS.ColumnCounter[ShieldingPAS.MatrixCounter];k++) … … 164 200 } 165 201 166 output << endl << "Total Times" << endl << "===============" << endl << Time.Header << endl; 167 for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) { 168 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) { 169 output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t"; 170 } 171 output << endl; 172 } 173 output << endl; 202 if (!NoTime) { 203 output << endl << "Total Times" << endl << "===============" << endl << Time.Header[Time.MatrixCounter] << endl; 204 for(int j=0;j<Time.RowCounter[Time.MatrixCounter];j++) { 205 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) { 206 output << scientific << Time.Matrix[ Time.MatrixCounter ][j][k] << "\t"; 207 } 208 output << endl; 209 } 210 output << endl; 211 } 174 212 output.close(); 175 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) 176 Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k]; 213 if (!NoTime) 214 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) 215 Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k] = Time.Matrix[ Time.MatrixCounter ][Time.RowCounter[Time.MatrixCounter]-1][k]; 177 216 178 217 // +++++++++++++++ ANALYZING ++++++++++++++++++++++++++++++ … … 184 223 // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order 185 224 // +++++++++++++++++++++++++++++++++++++++ Plotting Delta Simtime vs Bond Order 186 if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false; 187 if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false; 188 for(int j=Time.RowCounter[Time.MatrixCounter];j--;) 189 for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) { 190 Time.Matrix[ Time.MatrixCounter ][j][k] = 0.; 191 } 192 counter = 0; 193 output << "#Order\tFrag.No.\t" << Time.Header << endl; 194 output2 << "#Order\tFrag.No.\t" << Time.Header << endl; 195 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 196 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) 197 for(int j=Time.RowCounter[Time.MatrixCounter];j--;) 198 for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) { 199 Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k]; 200 } 201 counter += KeySet.FragmentsPerOrder[BondOrder]; 202 output << BondOrder+1 << "\t" << counter; 203 output2 << BondOrder+1 << "\t" << counter; 204 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) { 205 output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k]; 206 if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON) 207 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]; 208 else 209 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k]; 210 } 211 output << endl; 212 output2 << endl; 213 } 214 output.close(); 215 output2.close(); 225 if (!NoTime) { 226 if (!OpenOutputFile(output, argv[3], "SimTime-Order.dat" )) return false; 227 if (!OpenOutputFile(output2, argv[3], "DeltaSimTime-Order.dat" )) return false; 228 for(int j=Time.RowCounter[Time.MatrixCounter];j--;) 229 for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) { 230 Time.Matrix[ Time.MatrixCounter ][j][k] = 0.; 231 } 232 counter = 0; 233 output << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl; 234 output2 << "#Order\tFrag.No.\t" << Time.Header[Time.MatrixCounter] << endl; 235 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 236 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) 237 for(int j=Time.RowCounter[Time.MatrixCounter];j--;) 238 for(int k=Time.ColumnCounter[Time.MatrixCounter];k--;) { 239 Time.Matrix[ Time.MatrixCounter ][j][k] += Time.Matrix[ KeySet.OrderSet[BondOrder][i] ][j][k]; 240 } 241 counter += KeySet.FragmentsPerOrder[BondOrder]; 242 output << BondOrder+1 << "\t" << counter; 243 output2 << BondOrder+1 << "\t" << counter; 244 for(int k=0;k<Time.ColumnCounter[Time.MatrixCounter];k++) { 245 output << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k]; 246 if (fabs(Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]) > MYEPSILON) 247 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k] / Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter] ][k]; 248 else 249 output2 << "\t" << scientific << Time.Matrix[ Time.MatrixCounter ][ Time.RowCounter[Time.MatrixCounter]-1 ][k]; 250 } 251 output << endl; 252 output2 << endl; 253 } 254 output.close(); 255 output2.close(); 256 } 257 258 if (!NoHessian) { 259 // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in hessian to full QM 260 if (!CreateDataDeltaHessianOrderPerAtom(Hessian, HessianFragments, KeySet, argv[3], "DeltaHessian_xx-Order", "Plot of error between approximated hessian and full hessian versus the Bond Order", datum)) return 1; 261 262 // ++++++++++++++++++++++++++++++++++++++Plotting Hessian vs. Order 263 if (!CreateDataHessianOrderPerAtom(HessianFragments, KeySet, argv[3], "Hessian_xx-Order", "Plot of approximated hessian versus the Bond Order", datum)) return 1; 264 if (!AppendOutputFile(output, argv[3], "Hessian_xx-Order.dat" )) return false; 265 output << endl << "# Full" << endl; 266 for(int j=0;j<Hessian.RowCounter[Hessian.MatrixCounter];j++) { 267 output << j << "\t"; 268 for(int k=0;k<Hessian.ColumnCounter[Force.MatrixCounter];k++) 269 output << scientific << Hessian.Matrix[ Hessian.MatrixCounter ][j][k] << "\t"; 270 output << endl; 271 } 272 output.close(); 273 } 216 274 217 275 // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings … … 305 363 yrange.str("[1e-8:1e+1]"); 306 364 307 // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order 308 if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "", 1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1; 365 if (!NoTime) { 366 // +++++++++++++++++++++++++++++++++++++++ Plotting Simtime vs Bond Order 367 if (!CreatePlotOrder(Time, KeySet, argv[3], "SimTime-Order", 1, "below", "y", "", 1, 1, "bond order k", "Evaluation time [s]", Orderxrange.str().c_str(), "", "1" , "with linespoints", EnergyPlotLine)) return 1; 368 } 309 369 310 370 // +++++++++++++++++++++++++++++++++++++++ Plotting deviation in energy to full QM -
TabularUnified src/datacreator.cpp ¶
rf731ae reeec8f 63 63 cout << msg << endl; 64 64 output << "# " << msg << ", created on " << datum; 65 output << "#Order\tFrag.No.\t" << Fragments.Header << endl;65 output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 66 66 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 67 67 for(int i=KeySet.FragmentsPerOrder[BondOrder];i--;) { … … 96 96 cout << msg << endl; 97 97 output << "# " << msg << ", created on " << datum; 98 output << "#Order\tFrag.No.\t" << Fragments.Header << endl;98 output << "#Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 99 99 Fragments.SetLastMatrix(Energy.Matrix[Energy.MatrixCounter],0); 100 100 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { … … 133 133 cout << msg << endl; 134 134 output << "# " << msg << ", created on " << datum; 135 output << "# Order\tFrag.No.\t" << Fragments.Header << endl;135 output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 136 136 Fragments.SetLastMatrix(0.,0); 137 137 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { … … 165 165 cout << msg << endl; 166 166 output << "# " << msg << ", created on " << datum; 167 output << "# Order\tFrag.No.\t" << Fragments.Header << endl;167 output << "# Order\tFrag.No.\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 168 168 Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter],0); 169 169 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { … … 198 198 cout << msg << endl; 199 199 output << "# " << msg << ", created on " << datum; 200 output << "# AtomNo\t" << Fragments.Header << endl;200 output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 201 201 Fragments.SetLastMatrix(Force.Matrix[Force.MatrixCounter], 0); 202 202 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { … … 244 244 cout << msg << endl; 245 245 output << "# " << msg << ", created on " << datum; 246 output << "# AtomNo\t" << Fragments.Header << endl;246 output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 247 247 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 248 248 //cout << "Current order is " << BondOrder << "." << endl; … … 262 262 }; 263 263 264 265 /** Plot hessian error vs. vs atom vs. order. 266 * \param &Hessian HessianMatrix containing reference values (in MatrixCounter matrix) 267 * \param &Fragments HessianMatrix class containing matrix values 268 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order 269 * \param *prefix prefix in filename (without ending) 270 * \param *msg message to be place in first line as a comment 271 * \param *datum current date and time 272 * \return true if file was written successfully 273 */ 274 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 275 { 276 stringstream filename; 277 ofstream output; 278 double norm = 0.; 279 280 filename << prefix << ".dat"; 281 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 282 cout << msg << endl; 283 output << "# " << msg << ", created on " << datum; 284 output << "# AtomNo\t" << Fragments.Header[Fragments.MatrixCounter] << endl; 285 Fragments.SetLastMatrix(Hessian.Matrix[Hessian.MatrixCounter], 0); 286 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 287 //cout << "Current order is " << BondOrder << "." << endl; 288 Fragments.SumSubHessians(Fragments, KeySet, BondOrder, -1.); 289 // errors per atom 290 output << endl << "#Order\t" << BondOrder+1 << endl; 291 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 292 output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t"; 293 for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) { 294 output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t"; 295 } 296 output << endl; 297 } 298 output << endl; 299 } 300 output.close(); 301 return true; 302 }; 303 304 /** Plot hessian error vs. vs atom vs. order. 305 * \param &Fragments HessianMatrix class containing matrix values 306 * \param KeySet KeySetContainer class holding bond KeySetContainer::Order 307 * \param *prefix prefix in filename (without ending) 308 * \param *msg message to be place in first line as a comment 309 * \param *datum current date and time 310 * \return true if file was written successfully 311 */ 312 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum) 313 { 314 stringstream filename; 315 ofstream output; 316 317 filename << prefix << ".dat"; 318 if (!OpenOutputFile(output, dir, filename.str().c_str())) return false; 319 cout << msg << endl; 320 output << "# " << msg << ", created on " << datum; 321 output << "# AtomNo\t" << Fragments.Header[ Fragments.MatrixCounter ] << endl; 322 Fragments.SetLastMatrix(0., 0); 323 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 324 //cout << "Current order is " << BondOrder << "." << endl; 325 Fragments.SumSubHessians(Fragments, KeySet, BondOrder, 1.); 326 // errors per atom 327 output << endl << "#Order\t" << BondOrder+1 << endl; 328 for(int j=0;j<Fragments.RowCounter[ Fragments.MatrixCounter ];j++) { 329 output << Fragments.Indices[Fragments.MatrixCounter][j] << "\t"; 330 for (int l=0;l<Fragments.ColumnCounter[ Fragments.MatrixCounter ];l++) 331 output << scientific << Fragments.Matrix[Fragments.MatrixCounter][ j ][l] << "\t"; 332 output << endl; 333 } 334 output << endl; 335 } 336 output.close(); 337 return true; 338 }; 339 264 340 /** Plot matrix vs. fragment. 265 341 */ … … 273 349 cout << msg << endl; 274 350 output << "# " << msg << ", created on " << datum << endl; 275 output << "#Order\tFrag.No.\t" << Fragment.Header << endl;351 output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl; 276 352 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 277 353 for(int i=0;i<KeySet.FragmentsPerOrder[BondOrder];i++) { … … 338 414 cout << msg << endl; 339 415 output << "# " << msg << ", created on " << datum; 340 output << "#Order\tFrag.No.\t" << Fragment.Header << endl;416 output << "#Order\tFrag.No.\t" << Fragment.Header[ Fragment.MatrixCounter ] << endl; 341 417 // max 342 418 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { … … 538 614 void AbsEnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses) 539 615 { 540 stringstream line(Energy.Header );616 stringstream line(Energy.Header[ Energy.MatrixCounter ]); 541 617 string token; 542 618 … … 562 638 void EnergyPlotLine(ofstream &output, class MatrixContainer &Energy, const char *prefix, const char *xargument, const char *uses) 563 639 { 564 stringstream line(Energy.Header );640 stringstream line(Energy.Header[Energy.MatrixCounter]); 565 641 string token; 566 642 … … 586 662 void ForceMagnitudePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses) 587 663 { 588 stringstream line(Force.Header );664 stringstream line(Force.Header[Force.MatrixCounter]); 589 665 string token; 590 666 … … 617 693 void AbsFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses) 618 694 { 619 stringstream line(Force.Header );695 stringstream line(Force.Header[Force.MatrixCounter]); 620 696 string token; 621 697 … … 648 724 void BoxesForcePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses) 649 725 { 650 stringstream line(Force.Header );726 stringstream line(Force.Header[Force.MatrixCounter]); 651 727 char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"}; 652 728 string token; … … 680 756 void BoxesFirstForceValuePlotLine(ofstream &output, class MatrixContainer &Force, const char *prefix, const char *xargument, const char *uses) 681 757 { 682 stringstream line(Force.Header );758 stringstream line(Force.Header[Force.MatrixCounter]); 683 759 char *fillcolor[5] = {"black", "red", "blue", "green", "cyan"}; 684 760 string token; -
TabularUnified src/datacreator.hpp ¶
rf731ae reeec8f 26 26 bool CreateDataDeltaForcesOrder(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int)); 27 27 bool CreateDataDeltaForcesOrderPerAtom(class ForceMatrix &Force, class ForceMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 28 bool CreateDataHessianOrderPerAtom(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 29 bool CreateDataDeltaHessianOrderPerAtom(class HessianMatrix &Hessian, class HessianMatrix &Fragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum); 28 30 bool CreateDataFragment(class MatrixContainer &ForceFragments, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateForce)(class MatrixContainer &, int)); 29 31 bool CreateDataFragmentOrder(class MatrixContainer &Fragment, class KeySetsContainer &KeySet, char *dir, char *prefix, char *msg, char *datum, void (*CreateFragmentOrder)(class MatrixContainer &, class KeySetsContainer &, int)); -
TabularUnified src/joiner.cpp ¶
rf731ae reeec8f 37 37 stringstream prefix; 38 38 char *dir = NULL; 39 bool Hcorrected = true;39 bool NoHCorrection = false; 40 40 bool NoHessian = false; 41 41 … … 68 68 // ------------- Parse through all Fragment subdirs -------- 69 69 if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1; 70 Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0); 70 if (!Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0)) { 71 NoHCorrection = true; 72 cout << "No HCorrection matrices found, skipping these." << endl; 73 } 71 74 if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1; 72 75 if (!Hessian.ParseFragmentMatrix(argv[1], dir, HessianSuffix, 0,0)) { 73 76 NoHessian = true; 74 cout << "No hessian matrices found, skipping th ose.";77 cout << "No hessian matrices found, skipping these." << endl; 75 78 } 76 79 if (periode != NULL) { // also look for PAS values … … 81 84 // ---------- Parse the TE Factors into an array ----------------- 82 85 if (!Energy.InitialiseIndices()) return 1; 83 if (Hcorrected) Hcorrection.InitialiseIndices(); 86 if (!NoHCorrection) 87 Hcorrection.InitialiseIndices(); 84 88 85 89 // ---------- Parse the Force indices into an array --------------- … … 87 91 88 92 // ---------- Parse the Hessian (=force) indices into an array --------------- 89 if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1; 93 if (!NoHessian) 94 if (!Hessian.InitialiseIndices((class MatrixContainer *)&Force)) return 1; 90 95 91 96 // ---------- Parse the shielding indices into an array --------------- … … 100 105 if (!KeySet.ParseManyBodyTerms()) return 1; 101 106 if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1; 102 if (Hcorrected) HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter); 107 if (!NoHCorrection) 108 HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter); 103 109 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1; 104 110 if (!NoHessian) … … 126 132 cout << "Summing energy of order " << BondOrder+1 << " ..." << endl; 127 133 if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1; 128 if ( Hcorrected) {134 if (!NoHCorrection) { 129 135 HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder); 130 136 if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1; 131 if (Hcorrected)Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.);137 Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.); 132 138 } else 133 139 if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1; … … 159 165 if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1; 160 166 // hessian 161 if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1; 167 if (!NoHessian) 168 if (!Hessian.WriteLastMatrix(argv[1], (prefix.str()).c_str(), HessianSuffix)) return 1; 162 169 // shieldings 163 170 if (periode != NULL) { // also look for PAS values … … 170 177 prefix << dir << EnergyFragmentSuffix; 171 178 if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 172 if ( Hcorrected) {179 if (!NoHCorrection) { 173 180 prefix.str(" "); 174 181 prefix << dir << HcorrectionFragmentSuffix; … … 195 202 // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds 196 203 if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1; 197 if ( Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix);204 if (!NoHCorrection) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix); 198 205 if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1; 199 206 if (!NoHessian) -
TabularUnified src/parser.cpp ¶
rf731ae reeec8f 56 56 MatrixContainer::MatrixContainer() { 57 57 Indices = NULL; 58 Header = (char * ) Malloc(sizeof(char)*1023, "MatrixContainer::MatrixContainer:*Header");58 Header = (char **) Malloc(sizeof(char)*1, "MatrixContainer::MatrixContainer: **Header"); 59 59 Matrix = (double ***) Malloc(sizeof(double **)*(1), "MatrixContainer::MatrixContainer: ***Matrix"); // one more each for the total molecule 60 60 RowCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *RowCounter"); 61 61 ColumnCounter = (int *) Malloc(sizeof(int)*(1), "MatrixContainer::MatrixContainer: *ColumnCounter"); 62 Header[0] = NULL; 62 63 Matrix[0] = NULL; 63 64 RowCounter[0] = -1; … … 90 91 Free((void **)&Indices, "MatrixContainer::~MatrixContainer: **Indices"); 91 92 92 Free((void **)&Header, "MatrixContainer::~MatrixContainer: *Header"); 93 if (Header != NULL) 94 for(int i=MatrixCounter+1;i--;) 95 Free((void **)&Header[i], "MatrixContainer::~MatrixContainer: *Header[]"); 96 Free((void **)&Header, "MatrixContainer::~MatrixContainer: **Header"); 93 97 Free((void **)&RowCounter, "MatrixContainer::~MatrixContainer: *RowCounter"); 94 98 Free((void **)&ColumnCounter, "MatrixContainer::~MatrixContainer: *RowCounter"); … … 120 124 return false; 121 125 Indices[i] = (int *) Malloc(sizeof(int)*Matrix->RowCounter[i], "MatrixContainer::InitialiseIndices: *Indices[]"); 122 for(int j=Matrix->RowCounter[i];j--;) 126 for(int j=Matrix->RowCounter[i];j--;) { 123 127 Indices[i][j] = Matrix->Indices[i][j]; 128 //cout << Indices[i][j] << "\t"; 129 } 130 //cout << endl; 124 131 } 125 132 } … … 153 160 return false; 154 161 } 155 156 // skip some initial lines 162 163 // parse header 164 Header[MatrixNr] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseMatrix: *Header[]"); 157 165 for (int m=skiplines+1;m--;) 158 input.getline(Header , 1023);166 input.getline(Header[MatrixNr], 1023); 159 167 160 168 // scan header for number of columns 161 line.str(Header );169 line.str(Header[MatrixNr]); 162 170 for(int k=skipcolumns;k--;) 163 line >> Header ;171 line >> Header[MatrixNr]; 164 172 //cout << line.str() << endl; 165 173 ColumnCounter[MatrixNr]=0; … … 169 177 } 170 178 //cout << line.str() << endl; 171 cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;179 //cout << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl; 172 180 if (ColumnCounter[MatrixNr] == 0) 173 181 cerr << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl; … … 183 191 } 184 192 } 185 cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;193 //cout << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl; 186 194 if (RowCounter[MatrixNr] == 0) 187 195 cerr << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl; … … 189 197 // allocate matrix if it's not zero dimension in one direction 190 198 if ((ColumnCounter[MatrixNr] > 0) && (RowCounter[MatrixNr] > -1)) { 191 Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::Parse FragmentMatrix: **Matrix[]");199 Matrix[MatrixNr] = (double **) Malloc(sizeof(double *)*(RowCounter[MatrixNr]+1), "MatrixContainer::ParseMatrix: **Matrix[]"); 192 200 193 201 // parse in each entry for this matrix … … 195 203 input.seekg(ios::beg); 196 204 for (int m=skiplines+1;m--;) 197 input.getline(Header , 1023); // skip header198 line.str(Header );205 input.getline(Header[MatrixNr], 1023); // skip header 206 line.str(Header[MatrixNr]); 199 207 for(int k=skipcolumns;k--;) // skip columns in header too 200 208 line >> filename; 201 strncpy(Header , line.str().c_str(), 1023);209 strncpy(Header[MatrixNr], line.str().c_str(), 1023); 202 210 for(int j=0;j<RowCounter[MatrixNr];j++) { 203 Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::Parse FragmentMatrix: *Matrix[][]");211 Matrix[MatrixNr][j] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[][]"); 204 212 input.getline(filename, 1023); 205 213 stringstream lines(filename); … … 212 220 } 213 221 //cout << endl; 214 Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::Parse FragmentMatrix: *Matrix[RowCounter[MatrixNr]][]");222 Matrix[MatrixNr][ RowCounter[MatrixNr] ] = (double *) Malloc(sizeof(double)*ColumnCounter[MatrixNr], "MatrixContainer::ParseMatrix: *Matrix[RowCounter[MatrixNr]][]"); 215 223 for(int j=ColumnCounter[MatrixNr];j--;) 216 224 Matrix[MatrixNr][ RowCounter[MatrixNr] ][j] = 0.; … … 266 274 267 275 cout << "Parsing through each fragment and retrieving " << prefix << suffix << "." << endl; 276 Header = (char **) ReAlloc(Header, sizeof(char *)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: **Header"); // one more each for the total molecule 268 277 Matrix = (double ***) ReAlloc(Matrix, sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule 269 278 RowCounter = (int *) ReAlloc(RowCounter, sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter"); … … 271 280 for(int i=MatrixCounter+1;i--;) { 272 281 Matrix[i] = NULL; 282 Header[i] = NULL; 273 283 RowCounter[i] = -1; 274 284 ColumnCounter[i] = -1; … … 287 297 288 298 /** Allocates and resets the memory for a number \a MCounter of matrices. 289 * \param * GivenHeader Header line299 * \param **GivenHeader Header line for each matrix 290 300 * \param MCounter number of matrices 291 301 * \param *RCounter number of rows for each matrix 292 * \param *CCounter number of columns for each matri ces302 * \param *CCounter number of columns for each matrix 293 303 * \return Allocation successful 294 304 */ 295 bool MatrixContainer::AllocateMatrix(char *GivenHeader, int MCounter, int *RCounter, int *CCounter) 296 { 297 Header = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::ParseFragmentMatrix: *EnergyHeader"); 298 strncpy(Header, GivenHeader, 1023); 299 305 bool MatrixContainer::AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter) 306 { 300 307 MatrixCounter = MCounter; 301 Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: ***Matrix"); // one more each for the total molecule 302 RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *RowCounter"); 303 ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::ParseFragmentMatrix: *ColumnCounter"); 308 Header = (char **) Malloc(sizeof(char *)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *Header"); 309 Matrix = (double ***) Malloc(sizeof(double **)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: ***Matrix"); // one more each for the total molecule 310 RowCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *RowCounter"); 311 ColumnCounter = (int *) Malloc(sizeof(int)*(MatrixCounter+1), "MatrixContainer::AllocateMatrix: *ColumnCounter"); 304 312 for(int i=MatrixCounter+1;i--;) { 313 Header[i] = (char *) Malloc(sizeof(char)*1024, "MatrixContainer::AllocateMatrix: *Header[i]"); 314 strncpy(Header[i], GivenHeader[i], 1023); 305 315 RowCounter[i] = RCounter[i]; 306 316 ColumnCounter[i] = CCounter[i]; 307 Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer:: ParseFragmentMatrix: **Matrix[]");317 Matrix[i] = (double **) Malloc(sizeof(double *)*(RowCounter[i]+1), "MatrixContainer::AllocateMatrix: **Matrix[]"); 308 318 for(int j=RowCounter[i]+1;j--;) { 309 Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer:: ParseFragmentMatrix: *Matrix[][]");319 Matrix[i][j] = (double *) Malloc(sizeof(double)*ColumnCounter[i], "MatrixContainer::AllocateMatrix: *Matrix[][]"); 310 320 for(int k=ColumnCounter[i];k--;) 311 321 Matrix[i][j][k] = 0.; … … 385 395 }; 386 396 387 /** Sums the en ergy with each factor and put into last element of \a ***Energies.397 /** Sums the entries with each factor and put into last element of \a ***Matrix. 388 398 * Sums over "E"-terms to create the "F"-terms 389 399 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies. … … 462 472 return false; 463 473 } 464 output << Header << endl;474 output << Header[i] << endl; 465 475 for(int j=0;j<RowCounter[i];j++) { 466 476 for(int k=0;k<ColumnCounter[i];k++) … … 491 501 return false; 492 502 } 493 output << Header << endl;503 output << Header[MatrixCounter] << endl; 494 504 for(int j=0;j<RowCounter[MatrixCounter];j++) { 495 505 for(int k=0;k<ColumnCounter[MatrixCounter];k++) … … 765 775 int j = Indices[ FragmentNr ][l]; 766 776 if (j > RowCounter[MatrixCounter]) { 767 cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << " !" << endl;777 cerr << "Current hessian index " << j << " is greater than " << RowCounter[MatrixCounter] << ", where i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl; 768 778 return false; 769 779 } … … 772 782 int k = Indices[ FragmentNr ][m]; 773 783 if (k > ColumnCounter[MatrixCounter]) { 774 cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << " !" << endl;784 cerr << "Current hessian index " << k << " is greater than " << ColumnCounter[MatrixCounter] << ", where m=" << m << ", j=" << j << ", i=" << i << ", Order=" << Order << ", l=" << l << " and FragmentNr=" << FragmentNr << "!" << endl; 775 785 return false; 776 786 } … … 785 795 }; 786 796 797 /** Constructor for class HessianMatrix. 798 */ 799 HessianMatrix::HessianMatrix() : MatrixContainer() 800 { 801 IsSymmetric = true; 802 } 803 804 /** Sums the hessian entries with each factor and put into last element of \a ***Matrix. 805 * Sums over "E"-terms to create the "F"-terms 806 * \param Matrix MatrixContainer with matrices (LevelCounter by *ColumnCounter) with all the energies. 807 * \param KeySet KeySetContainer with bond Order and association mapping of each fragment to an order 808 * \param Order bond order 809 * \return true if summing was successful 810 */ 811 bool HessianMatrix::SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order) 812 { 813 // go through each order 814 for (int CurrentFragment=0;CurrentFragment<KeySet.FragmentsPerOrder[Order];CurrentFragment++) { 815 //cout << "Current Fragment is " << CurrentFragment << "/" << KeySet.OrderSet[Order][CurrentFragment] << "." << endl; 816 // then go per order through each suborder and pick together all the terms that contain this fragment 817 for(int SubOrder=0;SubOrder<=Order;SubOrder++) { // go through all suborders up to the desired order 818 for (int j=0;j<KeySet.FragmentsPerOrder[SubOrder];j++) { // go through all possible fragments of size suborder 819 if (KeySet.Contains(KeySet.OrderSet[Order][CurrentFragment], KeySet.OrderSet[SubOrder][j])) { 820 //cout << "Current other fragment is " << j << "/" << KeySet.OrderSet[SubOrder][j] << "." << endl; 821 // if the fragment's indices are all in the current fragment 822 for(int k=0;k<RowCounter[ KeySet.OrderSet[SubOrder][j] ];k++) { // go through all atoms in this fragment 823 int m = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][k]; 824 //cout << "Current row index is " << k << "/" << m << "." << endl; 825 if (m != -1) { // if it's not an added hydrogen 826 for (int l=0;l<RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ];l++) { // look for the corresponding index in the current fragment 827 //cout << "Comparing " << m << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l] << "." << endl; 828 if (m == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][l]) { 829 m = l; 830 break; 831 } 832 } 833 //cout << "Corresponding row index for " << k << " in CurrentFragment is " << m << "." << endl; 834 if (m > RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) { 835 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment] << " current row index " << m << " is greater than " << RowCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl; 836 return false; 837 } 838 839 for(int l=0;l<ColumnCounter[ KeySet.OrderSet[SubOrder][j] ];l++) { 840 int n = MatrixValues.Indices[ KeySet.OrderSet[SubOrder][j] ][l]; 841 //cout << "Current column index is " << l << "/" << n << "." << endl; 842 if (n != -1) { // if it's not an added hydrogen 843 for (int p=0;p<ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ];p++) { // look for the corresponding index in the current fragment 844 //cout << "Comparing " << n << " with " << MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p] << "." << endl; 845 if (n == MatrixValues.Indices[ KeySet.OrderSet[Order][CurrentFragment] ][p]) { 846 n = p; 847 break; 848 } 849 } 850 //cout << "Corresponding column index for " << l << " in CurrentFragment is " << n << "." << endl; 851 if (n > ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ]) { 852 cerr << "In fragment No. " << KeySet.OrderSet[Order][CurrentFragment] << " current column index " << n << " is greater than " << ColumnCounter[ KeySet.OrderSet[Order][CurrentFragment] ] << "!" << endl; 853 return false; 854 } 855 if (Order == SubOrder) { // equal order is always copy from Energies 856 Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] += MatrixValues.Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l]; 857 } else { 858 Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][m][n] -= Matrix[ KeySet.OrderSet[SubOrder][j] ][k][l]; 859 } 860 } 861 } 862 } 863 //if ((ColumnCounter[ KeySet.OrderSet[SubOrder][j] ]>1) && (RowCounter[0]-1 >= 1)) 864 //cout << "Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << RowCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][RowCounter[0]-1][1] << endl; 865 } 866 } else { 867 //cout << "Fragment " << KeySet.OrderSet[SubOrder][j] << " is not contained in fragment " << KeySet.OrderSet[Order][CurrentFragment] << "." << endl; 868 } 869 } 870 } 871 //cout << "Final Fragments[ KeySet.OrderSet[" << Order << "][" << CurrentFragment << "]=" << KeySet.OrderSet[Order][CurrentFragment] << " ][" << KeySet.AtomCounter[0]-1 << "][" << 1 << "] = " << Matrix[ KeySet.OrderSet[Order][CurrentFragment] ][KeySet.AtomCounter[0]-1][1] << endl; 872 } 873 874 return true; 875 }; 787 876 788 877 /** Calls MatrixContainer::ParseFragmentMatrix() and additionally allocates last plus one matrix. -
TabularUnified src/parser.hpp ¶
rf731ae reeec8f 45 45 double ***Matrix; 46 46 int **Indices; 47 char * Header;47 char **Header; 48 48 int MatrixCounter; 49 49 int *RowCounter; … … 56 56 bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr); 57 57 virtual bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns); 58 bool AllocateMatrix(char * GivenHeader, int MCounter, int *RCounter, int *CCounter);58 bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter); 59 59 bool ResetMatrix(); 60 60 double FindMinValue(); … … 88 88 // ======================================= CLASS HessianMatrix ============================= 89 89 90 class HessianMatrix : public MatrixContainer { 90 class HessianMatrix : public MatrixContainer { 91 91 public: 92 HessianMatrix(); 93 //~HessianMatrix(); 92 94 bool ParseIndices(char *name); 95 bool SumSubManyBodyTerms(class MatrixContainer &MatrixValues, class KeySetsContainer &KeySet, int Order); 93 96 bool SumSubHessians(class HessianMatrix &Fragments, class KeySetsContainer &KeySet, int Order, double sign); 94 97 bool ParseFragmentMatrix(char *name, char *prefix, string suffix, int skiplines, int skipcolumns); 98 private: 99 bool IsSymmetric; 95 100 }; 96 101
Note:
See TracChangeset
for help on using the changeset viewer.