Changeset d7e30c


Ignore:
Timestamp:
Aug 6, 2008, 8:59:06 AM (16 years ago)
Author:
Frederik Heber <heber@…>
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
Message:

Scaffold of new function: VerletForceIntegration() - does MD step by parsing external forces file.

PCP is too slow and unnecessarily ECut-dependent. If we incorporate the Verlet integration into molecuilder, we get a much more complete "wrapper" package.

Location:
src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • src/builder.cpp

    rb4b7c3 rd7e30c  
    147147          third->Output(1,2,(ofstream *)&cout);
    148148          fourth->Output(1,3,(ofstream *)&cout);
    149           n.MakeNormalVector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
    150           x.CopyVector(&second->x);
     149          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
     150          x.Copyvector(&second->x);
    151151          x.SubtractVector(&third->x);
    152           x.CopyVector(&fourth->x);
     152          x.Copyvector(&fourth->x);
    153153          x.SubtractVector(&third->x);
    154154         
     
    359359      } while ((param.type = periode->FindElement(shorthand)) == NULL);
    360360      cout << Verbose(0) << "Element is " << param.type->name << endl;
    361       mol->GetAlignVector(&param);
     361      mol->GetAlignvector(&param);
    362362      for (int i=NDIM;i--;) {
    363363        x.x[i] = gsl_vector_get(param.x,i);
     
    761761            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    762762            cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
    763             cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and 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;
    764764            cout << "\t-h/-H/-?\tGive this help screen." << endl;
    765765            cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
     
    767767            cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    768768            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;
    769770            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    770771            cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
     
    863864              if (!mol->AddXYZFile(argv[argptr]))
    864865                cout << Verbose(2) << "File not found." << endl;
    865               else
     866              else {
    866867                cout << Verbose(2) << "File found and parsed." << endl;
    867868                config_present = present;
     869              }
    868870              break;
    869871            default:   // no match? Don't step on (this is done in next switch's default)
     
    873875        if (config_present == present) {
    874876          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;
    875886            case 't':
    876887              ExitFlag = 1;
     
    959970              cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
    960971              break;
     972            case 'F':
    961973            case 'f':
    962974              if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
     
    10221034                int count;
    10231035                element ** Elements;
    1024                 vector ** Vectors;
     1036                vector ** vectors;
    10251037                if (faktor < 1) {
    10261038                  cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
     
    10311043                  count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
    10321044                  Elements = new element *[count];
    1033                   Vectors = new vector *[count];
     1045                  vectors = new vector *[count];
    10341046                  j = 0;
    10351047                  first = mol->start;
     
    10371049                    first = first->next;
    10381050                    Elements[j] = first->type;
    1039                     Vectors[j] = &first->x;
     1051                    vectors[j] = &first->x;
    10401052                    j++;
    10411053                  }
     
    10491061                    for (int k=count;k--;) { // go through every atom of the original cell
    10501062                      first = new atom(); // create a new body
    1051                       first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
     1063                      first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    10521064                      first->x.AddVector(&x);      // translate the coordinates
    10531065                      first->type = Elements[k];  // insert original element
     
    10571069                  // free memory
    10581070                  delete[](Elements);
    1059                   delete[](Vectors);
     1071                  delete[](vectors);
    10601072                  // correct cell size
    10611073                  if (axis < 0) { // if sign was negative, we have to translate everything
     
    11201132  clock_t start,end;
    11211133  element **Elements;
    1122   vector **Vectors;
     1134  vector **vectors;
    11231135
    11241136  // =========================== PARSE COMMAND LINE OPTIONS ====================================
     
    12191231          count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
    12201232          Elements = new element *[count];
    1221           Vectors = new vector *[count];
     1233          vectors = new vector *[count];
    12221234          j = 0;
    12231235          first = mol->start;
     
    12251237            first = first->next;
    12261238            Elements[j] = first->type;
    1227             Vectors[j] = &first->x;
     1239            vectors[j] = &first->x;
    12281240            j++;
    12291241          }
     
    12371249            for (int k=count;k--;) { // go through every atom of the original cell
    12381250              first = new atom(); // create a new body
    1239               first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
     1251              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
    12401252              first->x.AddVector(&x);      // translate the coordinates
    12411253              first->type = Elements[k];  // insert original element
     
    12471259          // free memory
    12481260          delete[](Elements);
    1249           delete[](Vectors);
     1261          delete[](vectors);
    12501262          // correct cell size
    12511263          if (axis < 0) { // if sign was negative, we have to translate everything
  • src/molecules.cpp

    rb4b7c3 rd7e30c  
    140140 * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
    141141 * -# 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::MakeNormalVector().
     142 *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
    143143 *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
    144144 *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
     
    151151 *    \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 )}}
    152152 *    \f]
    153  *    vector::GetNormalVector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalVector creates
     153 *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
    154154 *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
    155155 *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
     
    180180  atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
    181181  double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
    182   vector OrthoVector1, OrthoVector2;  // temporary vectors in coordination construction
    183   vector InBondVector;    // vector in direction of *Bond
     182  vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
     183  vector InBondvector;    // vector in direction of *Bond
    184184  bond *Binder = NULL;
    185185  double *matrix;
     
    187187//      *out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    188188  // create vector in direction of bond
    189   InBondVector.CopyVector(&TopReplacement->x);
    190   InBondVector.SubtractVector(&TopOrigin->x);
    191   bondlength = InBondVector.Norm();
     189  InBondvector.CopyVector(&TopReplacement->x);
     190  InBondvector.SubtractVector(&TopOrigin->x);
     191  bondlength = InBondvector.Norm();
    192192   
    193193   // 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 InBondVector have the wrong direction
     194   // the problem is not the H being out of the box, but InBondvector have the wrong direction
    195195   // due to TopReplacement or Origin being on the wrong side!
    196196  if (bondlength > BondDistance) {
    197 //    *out << Verbose(4) << "InBondVector is: ";
    198 //    InBondVector.Output(out);
     197//    *out << Verbose(4) << "InBondvector is: ";
     198//    InBondvector.Output(out);
    199199//    *out << endl;
    200     OrthoVector1.Zero();
     200    Orthovector1.Zero();
    201201    for (int i=NDIM;i--;) {
    202202      l = TopReplacement->x.x[i] - TopOrigin->x.x[i];
    203203      if (fabs(l) > BondDistance) { // is component greater than bond distance
    204         OrthoVector1.x[i] = (l < 0) ? -1. : +1.;
     204        Orthovector1.x[i] = (l < 0) ? -1. : +1.;
    205205      } // (signs are correct, was tested!)
    206206    }
    207207    matrix = ReturnFullMatrixforSymmetric(cell_size);
    208     OrthoVector1.MatrixMultiplication(matrix);
    209     InBondVector.SubtractVector(&OrthoVector1); // subtract just the additional translation
     208    Orthovector1.MatrixMultiplication(matrix);
     209    InBondvector.SubtractVector(&Orthovector1); // subtract just the additional translation
    210210    Free((void **)&matrix, "molecule::AddHydrogenReplacementAtom: *matrix");
    211     bondlength = InBondVector.Norm();
    212 //    *out << Verbose(4) << "Corrected InBondVector is now: ";
    213 //    InBondVector.Output(out);
     211    bondlength = InBondvector.Norm();
     212//    *out << Verbose(4) << "Corrected InBondvector is now: ";
     213//    InBondvector.Output(out);
    214214//    *out << endl;
    215215  } // periodic correction finished
    216216 
    217   InBondVector.Normalize();
     217  InBondvector.Normalize();
    218218  // get typical bond length and store as scale factor for later
    219219  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
     
    239239        FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
    240240      }
    241       InBondVector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
     241      InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
    242242      FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ...
    243       FirstOtherAtom->x.AddVector(&InBondVector);  // ... and add distance vector to replacement atom
     243      FirstOtherAtom->x.AddVector(&InBondvector);  // ... and add distance vector to replacement atom
    244244      AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
    245245//      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
     
    273273       
    274274        // determine the plane of these two with the *origin
    275         AllWentWell = AllWentWell && OrthoVector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
     275        AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
    276276      } else {
    277         OrthoVector1.GetOneNormalVector(&InBondVector);
     277        Orthovector1.GetOneNormalVector(&InBondvector);
    278278      }
    279279      //*out << Verbose(3)<< "Orthovector1: ";
    280       //OrthoVector1.Output(out);
     280      //Orthovector1.Output(out);
    281281      //*out << endl;
    282282      // orthogonal vector and bond vector between origin and replacement form the new plane
    283       OrthoVector1.MakeNormalVector(&InBondVector);
    284       OrthoVector1.Normalize();
    285       //*out << Verbose(3) << "ReScaleCheck: " << OrthoVector1.Norm() << " and " << InBondVector.Norm() << "." << endl;
     283      Orthovector1.MakeNormalVector(&InBondvector);
     284      Orthovector1.Normalize();
     285      //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
    286286     
    287287      // create the two Hydrogens ...
     
    302302      }
    303303      bondangle *= M_PI/180./2.;
    304 //      *out << Verbose(3) << "ReScaleCheck: InBondVector ";
    305 //      InBondVector.Output(out);
     304//      *out << Verbose(3) << "ReScaleCheck: InBondvector ";
     305//      InBondvector.Output(out);
    306306//      *out << endl;
    307307//      *out << Verbose(3) << "ReScaleCheck: Orthovector ";
    308 //      OrthoVector1.Output(out);
     308//      Orthovector1.Output(out);
    309309//      *out << endl;
    310310//      *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
    311311      FirstOtherAtom->x.Zero();
    312312      SecondOtherAtom->x.Zero();
    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));
     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));
    316316      }
    317317      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
     
    356356      ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
    357357
    358       // we need to vectors orthonormal the InBondVector
    359       AllWentWell = AllWentWell && OrthoVector1.GetOneNormalVector(&InBondVector);
     358      // we need to vectors orthonormal the InBondvector
     359      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector);
    360360//      *out << Verbose(3) << "Orthovector1: ";
    361 //      OrthoVector1.Output(out);
     361//      Orthovector1.Output(out);
    362362//      *out << endl;
    363       AllWentWell = AllWentWell && OrthoVector2.MakeNormalVector(&InBondVector, &OrthoVector1);
     363      AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
    364364//      *out << Verbose(3) << "Orthovector2: ";
    365 //      OrthoVector2.Output(out);
     365//      Orthovector2.Output(out);
    366366//      *out << endl;
    367367     
     
    371371      b = 2.*l*sin(alpha);    // base length of isosceles triangle
    372372      d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
    373       f = b/sqrt(3.);   // length for OrthVector1
    374       g = b/2.;         // length for OrthVector2
     373      f = b/sqrt(3.);   // length for Orthvector1
     374      g = b/2.;         // length for Orthvector2
    375375//      *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
    376376//      *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;
     
    378378      factors[1] = f;
    379379      factors[2] = 0.;
    380       FirstOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
     380      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    381381      factors[1] = -0.5*f;
    382382      factors[2] = g;
    383       SecondOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
     383      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    384384      factors[2] = -g;
    385       ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
     385      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
    386386
    387387      // rescale each to correct BondDistance
     
    802802  double tmp;
    803803  bool flag;
    804   vector TestVector, TranslationVector;
     804  vector Testvector, Translationvector;
    805805 
    806806  do {
     
    812812      if (Walker->type->Z != 1) {
    813813#endif
    814         TestVector.CopyVector(&Walker->x);
    815         TestVector.InverseMatrixMultiplication(matrix);
    816         TranslationVector.Zero();
     814        Testvector.CopyVector(&Walker->x);
     815        Testvector.InverseMatrixMultiplication(matrix);
     816        Translationvector.Zero();
    817817        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
    818818          Binder = ListOfBondsPerAtom[Walker->nr][i];
     
    824824                cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
    825825                if (tmp > 0)
    826                   TranslationVector.x[j] -= 1.;
     826                  Translationvector.x[j] -= 1.;
    827827                else
    828                   TranslationVector.x[j] += 1.;
     828                  Translationvector.x[j] += 1.;
    829829              }
    830830            }
    831831        }
    832         TestVector.AddVector(&TranslationVector);
    833         TestVector.MatrixMultiplication(matrix);
    834         Center.AddVector(&TestVector);
    835         cout << Verbose(1) << "Vector is: ";
    836         TestVector.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);
    837837        cout << endl;     
    838838#ifdef ADDHYDROGEN
     
    841841          Binder = ListOfBondsPerAtom[Walker->nr][i];
    842842          if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
    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);
     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);
    850850            cout << endl;     
    851851          }
     
    963963};
    964964
     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 */
     971bool 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
    9651059/** Align all atoms in such a manner that given vector \a *n is along z axis.
    9661060 * \param n[] alignment vector.
     
    11241218 * \bug this is not yet working properly it seems
    11251219 */
    1126 void molecule::GetAlignVector(struct lsq_params * par) const
     1220void molecule::GetAlignvector(struct lsq_params * par) const
    11271221{
    11281222    int np = 6;
     
    35733667          *out << (*runner) << " ";
    35743668        *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;
    35773671        InsertFragmentIntoGraph(out, FragmentSearch);
    35783672        //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
     
    38313925  enum Shading *ColorList = NULL;
    38323926  double tmp;
    3833   vector TranslationVector;
     3927  vector Translationvector;
    38343928  //class StackClass<atom *> *CompStack = NULL;
    38353929  class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
     
    38423936    // remove bonds that are beyond bonddistance
    38433937    for(int i=NDIM;i--;)
    3844       TranslationVector.x[i] = 0.;
     3938      Translationvector.x[i] = 0.;
    38453939    // scan all bonds
    38463940    Binder = first;
     
    38653959        tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
    38663960        if (fabs(tmp) > BondDistance)
    3867           TranslationVector.x[i] = (tmp < 0) ? +1. : -1.;
    3868       }
    3869       TranslationVector.MatrixMultiplication(matrix);
     3961          Translationvector.x[i] = (tmp < 0) ? +1. : -1.;
     3962      }
     3963      Translationvector.MatrixMultiplication(matrix);
    38703964      *out << Verbose(3) << "Translation vector is ";
    3871       TranslationVector.Output(out);
     3965      Translationvector.Output(out);
    38723966      *out << endl;
    38733967      // apply to all atoms of first component via BFS
     
    38793973        //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
    38803974        ColorList[Walker->nr] = black;    // mark as explored
    3881         Walker->x.AddVector(&TranslationVector); // translate
     3975        Walker->x.AddVector(&Translationvector); // translate
    38823976        for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
    38833977          if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
  • src/molecules.hpp

    rb4b7c3 rd7e30c  
    2020#include <set>
    2121#include <deque>
     22#include <list>
    2223
    2324#include "helpers.hpp"
     
    4647{
    4748  bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
     49};
     50
     51struct Trajectory
     52{
     53  double R[3];  //!< position vector
     54  double U[3];    //!< velocity vector
     55  double F[3];    //!< last force vector
    4856};
    4957
     
    219227        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
    220228        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
     229        bool VerletForceIntegration(char *file, double delta_t);
    221230       
    222231  bool CheckBounds(const vector *x) const;
    223   void GetAlignVector(struct lsq_params * par) const;
     232  void GetAlignvector(struct lsq_params * par) const;
    224233
    225234  /// Initialising routines in fragmentation 
Note: See TracChangeset for help on using the changeset viewer.