Changeset e08f45 for molecuilder/src/joiner.cpp
- Timestamp:
- Feb 9, 2009, 5:24:10 PM (17 years ago)
- Children:
- 451d7a
- Parents:
- 4aef8a
- git-author:
- Frederik Heber <heber@…> (02/09/09 15:55:37)
- git-committer:
- Frederik Heber <heber@…> (02/09/09 17:24:10)
- File:
-
- 1 edited
-
molecuilder/src/joiner.cpp (modified) (2 diffs, 1 prop)
Legend:
- Unmodified
- Added
- Removed
-
molecuilder/src/joiner.cpp
-
Property mode
changed from
100644to100755
r4aef8a re08f45 2 2 * 3 3 * Takes evaluated fragments (energy and forces) and by reading the factors files determines total energy 4 * and each force for the whole molecule. 5 * 4 * and each force for the whole molecule. 5 * 6 6 */ 7 7 8 8 //============================ INCLUDES =========================== 9 9 10 #include "datacreator.hpp" 11 #include "helpers.hpp" 12 #include "parser.hpp" 13 #include "periodentafel.hpp" 10 #include "datacreator.hpp" 11 #include "helpers.hpp" 12 #include "parser.hpp" 13 #include "periodentafel.hpp" 14 14 15 15 //============================== MAIN ============================= … … 17 17 int main(int argc, char **argv) 18 18 { 19 periodentafel *periode = NULL; // and a period table of all elements20 EnergyMatrix Energy;21 EnergyMatrix Hcorrection;22 ForceMatrix Force;23 EnergyMatrix EnergyFragments;24 EnergyMatrix HcorrectionFragments;25 ForceMatrix ForceFragments;26 ForceMatrix Shielding;27 ForceMatrix ShieldingPAS;28 ForceMatrix ShieldingFragments;29 ForceMatrix ShieldingPASFragments;30 KeySetsContainer KeySet; 31 stringstream prefix;32 char *dir = NULL;33 bool Hcorrected = true;19 periodentafel *periode = NULL; // and a period table of all elements 20 EnergyMatrix Energy; 21 EnergyMatrix Hcorrection; 22 ForceMatrix Force; 23 EnergyMatrix EnergyFragments; 24 EnergyMatrix HcorrectionFragments; 25 ForceMatrix ForceFragments; 26 ForceMatrix Shielding; 27 ForceMatrix ShieldingPAS; 28 ForceMatrix ShieldingFragments; 29 ForceMatrix ShieldingPASFragments; 30 KeySetsContainer KeySet; 31 stringstream prefix; 32 char *dir = NULL; 33 bool Hcorrected = true; 34 34 35 cout << "Joiner" << endl; 36 cout << "======" << endl; 37 38 // Get the command line options 39 if (argc < 3) { 40 cout << "Usage: " << argv[0] << " <inputdir> <prefix> [elementsdb]" << endl; 41 cout << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl; 42 cout << "<prefix>\tprefix of energy and forces file." << endl; 43 cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl; 44 return 1; 45 } else { 46 dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir"); 47 strcpy(dir, "/"); 48 strcat(dir, argv[2]); 49 } 50 if (argc > 3) { 51 periode = new periodentafel; 52 periode->LoadPeriodentafel(argv[3]); 53 } 54 55 // Test the given directory 56 if (!TestParams(argc, argv)) 57 return 1; 58 59 // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++ 60 61 // ------------- Parse through all Fragment subdirs -------- 62 if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1; 63 Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0); 64 if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1; 65 if (periode != NULL) { // also look for PAS values 66 if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1; 67 if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1; 68 } 35 cout << "Joiner" << endl; 36 cout << "======" << endl; 69 37 70 // ---------- Parse the TE Factors into an array ----------------- 71 if (!Energy.ParseIndices()) return 1; 72 if (Hcorrected) Hcorrection.ParseIndices(); 73 74 // ---------- Parse the Force indices into an array --------------- 75 if (!Force.ParseIndices(argv[1])) return 1; 38 // Get the command line options 39 if (argc < 3) { 40 cout << "Usage: " << argv[0] << " <inputdir> <prefix> [elementsdb]" << endl; 41 cout << "<inputdir>\ttherein the output of a molecuilder fragmentation is expected, each fragment with a subdir containing an energy.all and a forces.all file." << endl; 42 cout << "<prefix>\tprefix of energy and forces file." << endl; 43 cout << "[elementsdb]\tpath to elements database, needed for shieldings." << endl; 44 return 1; 45 } else { 46 dir = (char *) Malloc(sizeof(char)*(strlen(argv[2])+2), "main: *dir"); 47 strcpy(dir, "/"); 48 strcat(dir, argv[2]); 49 } 50 if (argc > 3) { 51 periode = new periodentafel; 52 periode->LoadPeriodentafel(argv[3]); 53 } 76 54 77 // ---------- Parse the shielding indices into an array --------------- 78 if (periode != NULL) { // also look for PAS values 79 if(!Shielding.ParseIndices(argv[1])) return 1; 80 if(!ShieldingPAS.ParseIndices(argv[1])) return 1; 81 } 55 // Test the given directory 56 if (!TestParams(argc, argv)) 57 return 1; 82 58 83 // ---------- Parse the KeySets into an array --------------- 84 if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1; 59 // +++++++++++++++++ PARSING +++++++++++++++++++++++++++++++ 85 60 86 if (!KeySet.ParseManyBodyTerms()) return 1; 87 if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1; 88 if (Hcorrected) HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter); 89 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1; 90 if (periode != NULL) { // also look for PAS values 91 if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1; 92 if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1; 93 } 94 95 // ----------- Resetting last matrices (where full QM values are stored right now) 96 if(!Energy.SetLastMatrix(0., 0)) return 1; 97 if(!Force.SetLastMatrix(0., 2)) return 1; 98 if (periode != NULL) { // also look for PAS values 99 if(!Shielding.SetLastMatrix(0., 2)) return 1; 100 if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1; 101 } 61 // ------------- Parse through all Fragment subdirs -------- 62 if (!Energy.ParseFragmentMatrix(argv[1], dir, EnergySuffix, 0,0)) return 1; 63 Hcorrected = Hcorrection.ParseFragmentMatrix(argv[1], "", HCORRECTIONSUFFIX, 0,0); 64 if (!Force.ParseFragmentMatrix(argv[1], dir, ForcesSuffix, 0,0)) return 1; 65 if (periode != NULL) { // also look for PAS values 66 if (!Shielding.ParseFragmentMatrix(argv[1], dir, ShieldingSuffix, 1, 0)) return 1; 67 if (!ShieldingPAS.ParseFragmentMatrix(argv[1], dir, ShieldingPASSuffix, 1, 0)) return 1; 68 } 102 69 103 // +++++++++++++++++ SUMMING +++++++++++++++++++++++++++++++ 70 // ---------- Parse the TE Factors into an array ----------------- 71 if (!Energy.ParseIndices()) return 1; 72 if (Hcorrected) Hcorrection.ParseIndices(); 104 73 105 // --------- sum up and write for each order---------------- 106 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 107 // --------- sum up energy -------------------- 108 cout << "Summing energy of order " << BondOrder+1 << " ..." << endl; 109 if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1; 110 if (Hcorrected) { 111 HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder); 112 if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1; 113 if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.); 114 } else 115 if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1; 116 // --------- sum up Forces -------------------- 117 cout << "Summing forces of order " << BondOrder+1 << " ..." << endl; 118 if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1; 119 if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1; 120 if (periode != NULL) { // also look for PAS values 121 cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl; 122 if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1; 123 if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1; 124 if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1; 125 if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1; 126 } 74 // ---------- Parse the Force indices into an array --------------- 75 if (!Force.ParseIndices(argv[1])) return 1; 127 76 128 // --------- write the energy and forces file -------------------- 129 prefix.str(" "); 130 prefix << dir << OrderSuffix << (BondOrder+1); 131 cout << "Writing files " << argv[1] << prefix.str() << ". ..." << endl; 132 // energy 133 if (!Energy.WriteLastMatrix(argv[1], (prefix.str()).c_str(), EnergySuffix)) return 1; 134 // forces 135 if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1; 136 // shieldings 137 if (periode != NULL) { // also look for PAS values 138 if (!Shielding.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingSuffix)) return 1; 139 if (!ShieldingPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingPASSuffix)) return 1; 140 } 141 } 142 // fragments 143 prefix.str(" "); 144 prefix << dir << EnergyFragmentSuffix; 145 if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 146 if (Hcorrected) { 147 prefix.str(" "); 148 prefix << dir << HcorrectionFragmentSuffix; 149 if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 150 } 151 prefix.str(" "); 152 prefix << dir << ForceFragmentSuffix; 153 if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 154 if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], FRAGMENTPREFIX ENERGYPERFRAGMENT, "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1; 155 if (periode != NULL) { // also look for PAS values 156 prefix.str(" "); 157 prefix << dir << ShieldingFragmentSuffix; 158 if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 159 prefix.str(" "); 160 prefix << dir << ShieldingPASFragmentSuffix; 161 if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 162 } 77 // ---------- Parse the shielding indices into an array --------------- 78 if (periode != NULL) { // also look for PAS values 79 if(!Shielding.ParseIndices(argv[1])) return 1; 80 if(!ShieldingPAS.ParseIndices(argv[1])) return 1; 81 } 163 82 164 // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds 165 if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1; 166 if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix); 167 if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1; 168 if (periode != NULL) { // also look for PAS values 169 if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1; 170 if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1; 171 } 83 // ---------- Parse the KeySets into an array --------------- 84 if (!KeySet.ParseKeySets(argv[1], Force.RowCounter, Force.MatrixCounter)) return 1; 172 85 173 // exit 174 delete(periode); 175 Free((void **)&dir, "main: *dir"); 176 cout << "done." << endl; 177 return 0; 86 if (!KeySet.ParseManyBodyTerms()) return 1; 87 if (!EnergyFragments.AllocateMatrix(Energy.Header, Energy.MatrixCounter, Energy.RowCounter, Energy.ColumnCounter)) return 1; 88 if (Hcorrected) HcorrectionFragments.AllocateMatrix(Hcorrection.Header, Hcorrection.MatrixCounter, Hcorrection.RowCounter, Hcorrection.ColumnCounter); 89 if (!ForceFragments.AllocateMatrix(Force.Header, Force.MatrixCounter, Force.RowCounter, Force.ColumnCounter)) return 1; 90 if (periode != NULL) { // also look for PAS values 91 if (!ShieldingFragments.AllocateMatrix(Shielding.Header, Shielding.MatrixCounter, Shielding.RowCounter, Shielding.ColumnCounter)) return 1; 92 if (!ShieldingPASFragments.AllocateMatrix(ShieldingPAS.Header, ShieldingPAS.MatrixCounter, ShieldingPAS.RowCounter, ShieldingPAS.ColumnCounter)) return 1; 93 } 94 95 // ----------- Resetting last matrices (where full QM values are stored right now) 96 if(!Energy.SetLastMatrix(0., 0)) return 1; 97 if(!Force.SetLastMatrix(0., 2)) return 1; 98 if (periode != NULL) { // also look for PAS values 99 if(!Shielding.SetLastMatrix(0., 2)) return 1; 100 if(!ShieldingPAS.SetLastMatrix(0., 2)) return 1; 101 } 102 103 // +++++++++++++++++ SUMMING +++++++++++++++++++++++++++++++ 104 105 // --------- sum up and write for each order---------------- 106 for (int BondOrder=0;BondOrder<KeySet.Order;BondOrder++) { 107 // --------- sum up energy -------------------- 108 cout << "Summing energy of order " << BondOrder+1 << " ..." << endl; 109 if (!EnergyFragments.SumSubManyBodyTerms(Energy, KeySet, BondOrder)) return 1; 110 if (Hcorrected) { 111 HcorrectionFragments.SumSubManyBodyTerms(Hcorrection, KeySet, BondOrder); 112 if (!Energy.SumSubEnergy(EnergyFragments, &HcorrectionFragments, KeySet, BondOrder, 1.)) return 1; 113 if (Hcorrected) Hcorrection.SumSubEnergy(HcorrectionFragments, NULL, KeySet, BondOrder, 1.); 114 } else 115 if (!Energy.SumSubEnergy(EnergyFragments, NULL, KeySet, BondOrder, 1.)) return 1; 116 // --------- sum up Forces -------------------- 117 cout << "Summing forces of order " << BondOrder+1 << " ..." << endl; 118 if (!ForceFragments.SumSubManyBodyTerms(Force, KeySet, BondOrder)) return 1; 119 if (!Force.SumSubForces(ForceFragments, KeySet, BondOrder, 1.)) return 1; 120 if (periode != NULL) { // also look for PAS values 121 cout << "Summing shieldings of order " << BondOrder+1 << " ..." << endl; 122 if (!ShieldingFragments.SumSubManyBodyTerms(Shielding, KeySet, BondOrder)) return 1; 123 if (!Shielding.SumSubForces(ShieldingFragments, KeySet, BondOrder, 1.)) return 1; 124 if (!ShieldingPASFragments.SumSubManyBodyTerms(ShieldingPAS, KeySet, BondOrder)) return 1; 125 if (!ShieldingPAS.SumSubForces(ShieldingPASFragments, KeySet, BondOrder, 1.)) return 1; 126 } 127 128 // --------- write the energy and forces file -------------------- 129 prefix.str(" "); 130 prefix << dir << OrderSuffix << (BondOrder+1); 131 cout << "Writing files " << argv[1] << prefix.str() << ". ..." << endl; 132 // energy 133 if (!Energy.WriteLastMatrix(argv[1], (prefix.str()).c_str(), EnergySuffix)) return 1; 134 // forces 135 if (!Force.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ForcesSuffix)) return 1; 136 // shieldings 137 if (periode != NULL) { // also look for PAS values 138 if (!Shielding.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingSuffix)) return 1; 139 if (!ShieldingPAS.WriteLastMatrix(argv[1], (prefix.str()).c_str(), ShieldingPASSuffix)) return 1; 140 } 141 } 142 // fragments 143 prefix.str(" "); 144 prefix << dir << EnergyFragmentSuffix; 145 if (!EnergyFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 146 if (Hcorrected) { 147 prefix.str(" "); 148 prefix << dir << HcorrectionFragmentSuffix; 149 if (!HcorrectionFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 150 } 151 prefix.str(" "); 152 prefix << dir << ForceFragmentSuffix; 153 if (!ForceFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 154 if (!CreateDataFragment(EnergyFragments, KeySet, argv[1], FRAGMENTPREFIX ENERGYPERFRAGMENT, "fragment energy versus the Fragment No", "today", CreateEnergy)) return 1; 155 if (periode != NULL) { // also look for PAS values 156 prefix.str(" "); 157 prefix << dir << ShieldingFragmentSuffix; 158 if (!ShieldingFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 159 prefix.str(" "); 160 prefix << dir << ShieldingPASFragmentSuffix; 161 if (!ShieldingPASFragments.WriteTotalFragments(argv[1], (prefix.str()).c_str())) return 1; 162 } 163 164 // write last matrices as fragments into central dir (not subdir as above), for analyzer to know index bounds 165 if (!Energy.WriteLastMatrix(argv[1], dir, EnergyFragmentSuffix)) return 1; 166 if (Hcorrected) Hcorrection.WriteLastMatrix(argv[1], dir, HcorrectionFragmentSuffix); 167 if (!Force.WriteLastMatrix(argv[1], dir, ForceFragmentSuffix)) return 1; 168 if (periode != NULL) { // also look for PAS values 169 if (!Shielding.WriteLastMatrix(argv[1], dir, ShieldingFragmentSuffix)) return 1; 170 if (!ShieldingPAS.WriteLastMatrix(argv[1], dir, ShieldingPASFragmentSuffix)) return 1; 171 } 172 173 // exit 174 delete(periode); 175 Free((void **)&dir, "main: *dir"); 176 cout << "done." << endl; 177 return 0; 178 178 }; 179 179 -
Property mode
changed from
Note:
See TracChangeset
for help on using the changeset viewer.
