Changeset d7e30c
- Timestamp:
- Aug 6, 2008, 8:59:06 AM (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:
- e9b8bb
- Parents:
- b4b7c3
- Location:
- src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
rb4b7c3 rd7e30c 147 147 third->Output(1,2,(ofstream *)&cout); 148 148 fourth->Output(1,3,(ofstream *)&cout); 149 n.MakeNormal Vector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);150 x.Copy Vector(&second->x);149 n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x); 150 x.Copyvector(&second->x); 151 151 x.SubtractVector(&third->x); 152 x.Copy Vector(&fourth->x);152 x.Copyvector(&fourth->x); 153 153 x.SubtractVector(&third->x); 154 154 … … 359 359 } while ((param.type = periode->FindElement(shorthand)) == NULL); 360 360 cout << Verbose(0) << "Element is " << param.type->name << endl; 361 mol->GetAlign Vector(¶m);361 mol->GetAlignvector(¶m); 362 362 for (int i=NDIM;i--;) { 363 363 x.x[i] = gsl_vector_get(param.x,i); … … 761 761 cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl; 762 762 cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl; 763 cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA mannerand stores config files in same dir as config." << endl;763 cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl; 764 764 cout << "\t-h/-H/-?\tGive this help screen." << endl; 765 765 cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl; … … 767 767 cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl; 768 768 cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl; 769 cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl; 769 770 cout << "\t-r\t\tConvert file from an old pcp syntax." << endl; 770 771 cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl; … … 863 864 if (!mol->AddXYZFile(argv[argptr])) 864 865 cout << Verbose(2) << "File not found." << endl; 865 else 866 else { 866 867 cout << Verbose(2) << "File found and parsed." << endl; 867 868 config_present = present; 869 } 868 870 break; 869 871 default: // no match? Don't step on (this is done in next switch's default) … … 873 875 if (config_present == present) { 874 876 switch(argv[argptr-1][1]) { 877 case 'P': 878 ExitFlag = 1; 879 cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl; 880 if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1]))) 881 cout << Verbose(2) << "File not found." << endl; 882 else 883 cout << Verbose(2) << "File found and parsed." << endl; 884 argptr+=2; 885 break; 875 886 case 't': 876 887 ExitFlag = 1; … … 959 970 cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl; 960 971 break; 972 case 'F': 961 973 case 'f': 962 974 if (ExitFlag ==0) ExitFlag = 2; // only set if not already by other command line switch … … 1022 1034 int count; 1023 1035 element ** Elements; 1024 vector ** Vectors;1036 vector ** vectors; 1025 1037 if (faktor < 1) { 1026 1038 cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl; … … 1031 1043 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand 1032 1044 Elements = new element *[count]; 1033 Vectors = new vector *[count];1045 vectors = new vector *[count]; 1034 1046 j = 0; 1035 1047 first = mol->start; … … 1037 1049 first = first->next; 1038 1050 Elements[j] = first->type; 1039 Vectors[j] = &first->x;1051 vectors[j] = &first->x; 1040 1052 j++; 1041 1053 } … … 1049 1061 for (int k=count;k--;) { // go through every atom of the original cell 1050 1062 first = new atom(); // create a new body 1051 first->x.CopyVector( Vectors[k]); // use coordinate of original atom1063 first->x.CopyVector(vectors[k]); // use coordinate of original atom 1052 1064 first->x.AddVector(&x); // translate the coordinates 1053 1065 first->type = Elements[k]; // insert original element … … 1057 1069 // free memory 1058 1070 delete[](Elements); 1059 delete[]( Vectors);1071 delete[](vectors); 1060 1072 // correct cell size 1061 1073 if (axis < 0) { // if sign was negative, we have to translate everything … … 1120 1132 clock_t start,end; 1121 1133 element **Elements; 1122 vector ** Vectors;1134 vector **vectors; 1123 1135 1124 1136 // =========================== PARSE COMMAND LINE OPTIONS ==================================== … … 1219 1231 count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand 1220 1232 Elements = new element *[count]; 1221 Vectors = new vector *[count];1233 vectors = new vector *[count]; 1222 1234 j = 0; 1223 1235 first = mol->start; … … 1225 1237 first = first->next; 1226 1238 Elements[j] = first->type; 1227 Vectors[j] = &first->x;1239 vectors[j] = &first->x; 1228 1240 j++; 1229 1241 } … … 1237 1249 for (int k=count;k--;) { // go through every atom of the original cell 1238 1250 first = new atom(); // create a new body 1239 first->x.CopyVector( Vectors[k]); // use coordinate of original atom1251 first->x.CopyVector(vectors[k]); // use coordinate of original atom 1240 1252 first->x.AddVector(&x); // translate the coordinates 1241 1253 first->type = Elements[k]; // insert original element … … 1247 1259 // free memory 1248 1260 delete[](Elements); 1249 delete[]( Vectors);1261 delete[](vectors); 1250 1262 // correct cell size 1251 1263 if (axis < 0) { // if sign was negative, we have to translate everything -
src/molecules.cpp
rb4b7c3 rd7e30c 140 140 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one 141 141 * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of 142 * *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormal Vector().142 * *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector(). 143 143 * The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two 144 144 * replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the … … 151 151 * \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}} 152 152 * \f] 153 * vector::GetNormal Vector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalVector creates153 * vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates 154 154 * the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above. 155 155 * The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that … … 180 180 atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added 181 181 double b,l,d,f,g, alpha, factors[NDIM]; // hold temporary values in triple bond case for coordination determination 182 vector Ortho Vector1, OrthoVector2; // temporary vectors in coordination construction183 vector InBond Vector; // vector in direction of *Bond182 vector Orthovector1, Orthovector2; // temporary vectors in coordination construction 183 vector InBondvector; // vector in direction of *Bond 184 184 bond *Binder = NULL; 185 185 double *matrix; … … 187 187 // *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl; 188 188 // create vector in direction of bond 189 InBond Vector.CopyVector(&TopReplacement->x);190 InBond Vector.SubtractVector(&TopOrigin->x);191 bondlength = InBond Vector.Norm();189 InBondvector.CopyVector(&TopReplacement->x); 190 InBondvector.SubtractVector(&TopOrigin->x); 191 bondlength = InBondvector.Norm(); 192 192 193 193 // is greater than typical bond distance? Then we have to correct periodically 194 // the problem is not the H being out of the box, but InBond Vector have the wrong direction194 // the problem is not the H being out of the box, but InBondvector have the wrong direction 195 195 // due to TopReplacement or Origin being on the wrong side! 196 196 if (bondlength > BondDistance) { 197 // *out << Verbose(4) << "InBond Vector is: ";198 // InBond Vector.Output(out);197 // *out << Verbose(4) << "InBondvector is: "; 198 // InBondvector.Output(out); 199 199 // *out << endl; 200 Ortho Vector1.Zero();200 Orthovector1.Zero(); 201 201 for (int i=NDIM;i--;) { 202 202 l = TopReplacement->x.x[i] - TopOrigin->x.x[i]; 203 203 if (fabs(l) > BondDistance) { // is component greater than bond distance 204 Ortho Vector1.x[i] = (l < 0) ? -1. : +1.;204 Orthovector1.x[i] = (l < 0) ? -1. : +1.; 205 205 } // (signs are correct, was tested!) 206 206 } 207 207 matrix = ReturnFullMatrixforSymmetric(cell_size); 208 Ortho Vector1.MatrixMultiplication(matrix);209 InBond Vector.SubtractVector(&OrthoVector1); // subtract just the additional translation208 Orthovector1.MatrixMultiplication(matrix); 209 InBondvector.SubtractVector(&Orthovector1); // subtract just the additional translation 210 210 Free((void **)&matrix, "molecule::AddHydrogenReplacementAtom: *matrix"); 211 bondlength = InBond Vector.Norm();212 // *out << Verbose(4) << "Corrected InBond Vector is now: ";213 // InBond Vector.Output(out);211 bondlength = InBondvector.Norm(); 212 // *out << Verbose(4) << "Corrected InBondvector is now: "; 213 // InBondvector.Output(out); 214 214 // *out << endl; 215 215 } // periodic correction finished 216 216 217 InBond Vector.Normalize();217 InBondvector.Normalize(); 218 218 // get typical bond length and store as scale factor for later 219 219 BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1]; … … 239 239 FirstOtherAtom->father = NULL; // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father 240 240 } 241 InBond Vector.Scale(&BondRescale); // rescale the distance vector to Hydrogen bond length241 InBondvector.Scale(&BondRescale); // rescale the distance vector to Hydrogen bond length 242 242 FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ... 243 FirstOtherAtom->x.AddVector(&InBond Vector); // ... and add distance vector to replacement atom243 FirstOtherAtom->x.AddVector(&InBondvector); // ... and add distance vector to replacement atom 244 244 AllWentWell = AllWentWell && AddAtom(FirstOtherAtom); 245 245 // *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: "; … … 273 273 274 274 // determine the plane of these two with the *origin 275 AllWentWell = AllWentWell && Ortho Vector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);275 AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x); 276 276 } else { 277 Ortho Vector1.GetOneNormalVector(&InBondVector);277 Orthovector1.GetOneNormalVector(&InBondvector); 278 278 } 279 279 //*out << Verbose(3)<< "Orthovector1: "; 280 //Ortho Vector1.Output(out);280 //Orthovector1.Output(out); 281 281 //*out << endl; 282 282 // orthogonal vector and bond vector between origin and replacement form the new plane 283 Ortho Vector1.MakeNormalVector(&InBondVector);284 Ortho Vector1.Normalize();285 //*out << Verbose(3) << "ReScaleCheck: " << Ortho Vector1.Norm() << " and " << InBondVector.Norm() << "." << endl;283 Orthovector1.MakeNormalVector(&InBondvector); 284 Orthovector1.Normalize(); 285 //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl; 286 286 287 287 // create the two Hydrogens ... … … 302 302 } 303 303 bondangle *= M_PI/180./2.; 304 // *out << Verbose(3) << "ReScaleCheck: InBond Vector ";305 // InBond Vector.Output(out);304 // *out << Verbose(3) << "ReScaleCheck: InBondvector "; 305 // InBondvector.Output(out); 306 306 // *out << endl; 307 307 // *out << Verbose(3) << "ReScaleCheck: Orthovector "; 308 // Ortho Vector1.Output(out);308 // Orthovector1.Output(out); 309 309 // *out << endl; 310 310 // *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl; 311 311 FirstOtherAtom->x.Zero(); 312 312 SecondOtherAtom->x.Zero(); 313 for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBond Vector is bondangle = 0 direction)314 FirstOtherAtom->x.x[i] = InBond Vector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (sin(bondangle));315 SecondOtherAtom->x.x[i] = InBond Vector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (-sin(bondangle));313 for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction) 314 FirstOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (sin(bondangle)); 315 SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle)); 316 316 } 317 317 FirstOtherAtom->x.Scale(&BondRescale); // rescale by correct BondDistance … … 356 356 ThirdOtherAtom->father = NULL; // we are just an added hydrogen with no father 357 357 358 // we need to vectors orthonormal the InBond Vector359 AllWentWell = AllWentWell && Ortho Vector1.GetOneNormalVector(&InBondVector);358 // we need to vectors orthonormal the InBondvector 359 AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector); 360 360 // *out << Verbose(3) << "Orthovector1: "; 361 // Ortho Vector1.Output(out);361 // Orthovector1.Output(out); 362 362 // *out << endl; 363 AllWentWell = AllWentWell && Ortho Vector2.MakeNormalVector(&InBondVector, &OrthoVector1);363 AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1); 364 364 // *out << Verbose(3) << "Orthovector2: "; 365 // Ortho Vector2.Output(out);365 // Orthovector2.Output(out); 366 366 // *out << endl; 367 367 … … 371 371 b = 2.*l*sin(alpha); // base length of isosceles triangle 372 372 d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.); // length for InBondvector 373 f = b/sqrt(3.); // length for Orth Vector1374 g = b/2.; // length for Orth Vector2373 f = b/sqrt(3.); // length for Orthvector1 374 g = b/2.; // length for Orthvector2 375 375 // *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl; 376 376 // *out << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl; … … 378 378 factors[1] = f; 379 379 factors[2] = 0.; 380 FirstOtherAtom->x.LinearCombinationOfVectors(&InBond Vector, &OrthoVector1, &OrthoVector2, factors);380 FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors); 381 381 factors[1] = -0.5*f; 382 382 factors[2] = g; 383 SecondOtherAtom->x.LinearCombinationOfVectors(&InBond Vector, &OrthoVector1, &OrthoVector2, factors);383 SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors); 384 384 factors[2] = -g; 385 ThirdOtherAtom->x.LinearCombinationOfVectors(&InBond Vector, &OrthoVector1, &OrthoVector2, factors);385 ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors); 386 386 387 387 // rescale each to correct BondDistance … … 802 802 double tmp; 803 803 bool flag; 804 vector Test Vector, TranslationVector;804 vector Testvector, Translationvector; 805 805 806 806 do { … … 812 812 if (Walker->type->Z != 1) { 813 813 #endif 814 Test Vector.CopyVector(&Walker->x);815 Test Vector.InverseMatrixMultiplication(matrix);816 Translation Vector.Zero();814 Testvector.CopyVector(&Walker->x); 815 Testvector.InverseMatrixMultiplication(matrix); 816 Translationvector.Zero(); 817 817 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) { 818 818 Binder = ListOfBondsPerAtom[Walker->nr][i]; … … 824 824 cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl; 825 825 if (tmp > 0) 826 Translation Vector.x[j] -= 1.;826 Translationvector.x[j] -= 1.; 827 827 else 828 Translation Vector.x[j] += 1.;828 Translationvector.x[j] += 1.; 829 829 } 830 830 } 831 831 } 832 Test Vector.AddVector(&TranslationVector);833 Test Vector.MatrixMultiplication(matrix);834 Center.AddVector(&Test Vector);835 cout << Verbose(1) << " Vector is: ";836 Test Vector.Output((ofstream *)&cout);832 Testvector.AddVector(&Translationvector); 833 Testvector.MatrixMultiplication(matrix); 834 Center.AddVector(&Testvector); 835 cout << Verbose(1) << "vector is: "; 836 Testvector.Output((ofstream *)&cout); 837 837 cout << endl; 838 838 #ifdef ADDHYDROGEN … … 841 841 Binder = ListOfBondsPerAtom[Walker->nr][i]; 842 842 if (Binder->GetOtherAtom(Walker)->type->Z == 1) { 843 Test Vector.CopyVector(&Binder->GetOtherAtom(Walker)->x);844 Test Vector.InverseMatrixMultiplication(matrix);845 Test Vector.AddVector(&TranslationVector);846 Test Vector.MatrixMultiplication(matrix);847 Center.AddVector(&Test Vector);848 cout << Verbose(1) << "Hydrogen Vector is: ";849 Test Vector.Output((ofstream *)&cout);843 Testvector.CopyVector(&Binder->GetOtherAtom(Walker)->x); 844 Testvector.InverseMatrixMultiplication(matrix); 845 Testvector.AddVector(&Translationvector); 846 Testvector.MatrixMultiplication(matrix); 847 Center.AddVector(&Testvector); 848 cout << Verbose(1) << "Hydrogen vector is: "; 849 Testvector.Output((ofstream *)&cout); 850 850 cout << endl; 851 851 } … … 963 963 }; 964 964 965 /** Parses nuclear forces from file and performs Verlet integration. 966 * This adds a new MD step to the config file. 967 * \param *file filename 968 * \param delta_t time step width in atomic units 969 * \return true - file found and parsed, false - file not found or imparsable 970 */ 971 bool molecule::VerletForceIntegration(char *file, double delta_t) 972 { 973 bool status = true; 974 element *runner = elemente->start; 975 atom *walker = NULL; 976 int ElementNo, AtomNo; 977 ifstream input(file); 978 double Forcesvector[3]; 979 double dummy; 980 string token; 981 size_t oldpos, newpos; 982 stringstream item; 983 int i, Ion[2]; 984 double a, IonMass; 985 986 CountElements(); // make sure ElementsInMolecule is up to date 987 988 // check file 989 if (input == NULL) { 990 return false; 991 } else { 992 ElementNo = 0; 993 while (runner->next != elemente->end) { // go through every element 994 runner = runner->next; 995 IonMass = runner->mass; 996 if (ElementsInMolecule[runner->Z]) { // if this element got atoms 997 ElementNo++; 998 AtomNo = 0; 999 walker = start; 1000 while (walker->next != end) { // go through every atom of this element 1001 walker = walker->next; 1002 if (walker->type == runner) { // if this atom fits to element 1003 AtomNo++; 1004 // parse Forcevector for this ion 1005 getline (input, token, '\n'); 1006 newpos = oldpos = i = 0; 1007 while ((newpos = token.find( '\t', oldpos)) != oldpos) { 1008 i++; 1009 item.str(); 1010 switch (i) { 1011 case 1: 1012 case 2: 1013 item >> Ion[i-1]; 1014 break; 1015 case 6: 1016 case 7: 1017 case 8: 1018 item >> Forcesvector[i-6]; 1019 break; 1020 case 3: 1021 if ((Ion[0] != ElementNo) || (Ion[1] != AtomNo)) 1022 cout << "ERROR: Expected " << ElementNo << ", " << AtomNo << "but parsed " << Ion[0] << ", " << Ion[1] << "." << endl; 1023 // still parse into dummy! Hence no break here ... 1024 default: 1025 item >> dummy; 1026 } 1027 oldpos = newpos; 1028 } 1029 1030 // perform Verlet integration for this atom with position, velocity and force vector 1031 1032 // UpdateR 1033 // for (int d=0; d<NDIM; d++) { 1034 // R_old_old[d] = R_old[d]; // shift old values 1035 // R_old[d] = R[d]; 1036 // R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2 1037 // FIon_old[d] = FIon[d]; // store old force 1038 // } 1039 // // UpdateU 1040 // a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a 1041 // for (int d=0; d<NDIM; d++) { 1042 // U[d] += a * (FIon[d]+FIon_old[d]); 1043 // } 1044 1045 1046 // add trajectory point 1047 1048 } 1049 } 1050 } 1051 } 1052 return true; 1053 } 1054 1055 // exit 1056 return status; 1057 }; 1058 965 1059 /** Align all atoms in such a manner that given vector \a *n is along z axis. 966 1060 * \param n[] alignment vector. … … 1124 1218 * \bug this is not yet working properly it seems 1125 1219 */ 1126 void molecule::GetAlign Vector(struct lsq_params * par) const1220 void molecule::GetAlignvector(struct lsq_params * par) const 1127 1221 { 1128 1222 int np = 6; … … 3573 3667 *out << (*runner) << " "; 3574 3668 *out << endl; 3575 if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))3576 *out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;3669 //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet)) 3670 //*out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl; 3577 3671 InsertFragmentIntoGraph(out, FragmentSearch); 3578 3672 //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList); … … 3831 3925 enum Shading *ColorList = NULL; 3832 3926 double tmp; 3833 vector Translation Vector;3927 vector Translationvector; 3834 3928 //class StackClass<atom *> *CompStack = NULL; 3835 3929 class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount); … … 3842 3936 // remove bonds that are beyond bonddistance 3843 3937 for(int i=NDIM;i--;) 3844 Translation Vector.x[i] = 0.;3938 Translationvector.x[i] = 0.; 3845 3939 // scan all bonds 3846 3940 Binder = first; … … 3865 3959 tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]; 3866 3960 if (fabs(tmp) > BondDistance) 3867 Translation Vector.x[i] = (tmp < 0) ? +1. : -1.;3868 } 3869 Translation Vector.MatrixMultiplication(matrix);3961 Translationvector.x[i] = (tmp < 0) ? +1. : -1.; 3962 } 3963 Translationvector.MatrixMultiplication(matrix); 3870 3964 *out << Verbose(3) << "Translation vector is "; 3871 Translation Vector.Output(out);3965 Translationvector.Output(out); 3872 3966 *out << endl; 3873 3967 // apply to all atoms of first component via BFS … … 3879 3973 //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl; 3880 3974 ColorList[Walker->nr] = black; // mark as explored 3881 Walker->x.AddVector(&Translation Vector); // translate3975 Walker->x.AddVector(&Translationvector); // translate 3882 3976 for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) { // go through all binding partners 3883 3977 if (ListOfBondsPerAtom[Walker->nr][i] != Binder) { -
src/molecules.hpp
rb4b7c3 rd7e30c 20 20 #include <set> 21 21 #include <deque> 22 #include <list> 22 23 23 24 #include "helpers.hpp" … … 46 47 { 47 48 bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const; 49 }; 50 51 struct Trajectory 52 { 53 double R[3]; //!< position vector 54 double U[3]; //!< velocity vector 55 double F[3]; //!< last force vector 48 56 }; 49 57 … … 219 227 void PrincipalAxisSystem(ofstream *out, bool DoRotate); 220 228 double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem); 229 bool VerletForceIntegration(char *file, double delta_t); 221 230 222 231 bool CheckBounds(const vector *x) const; 223 void GetAlign Vector(struct lsq_params * par) const;232 void GetAlignvector(struct lsq_params * par) const; 224 233 225 234 /// Initialising routines in fragmentation
Note:
See TracChangeset
for help on using the changeset viewer.