Ignore:
Timestamp:
Apr 29, 2010, 1:55:21 PM (16 years ago)
Author:
Tillmann Crueger <crueger@…>
Children:
070651, 5d1a94
Parents:
90c4460 (diff), 32842d8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'VectorRefactoring' into StructureRefactoring

Conflicts:

molecuilder/src/Legacy/oldmenu.cpp
molecuilder/src/Makefile.am
molecuilder/src/analysis_correlation.cpp
molecuilder/src/boundary.cpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/molecule.cpp
molecuilder/src/molecule_fragmentation.cpp
molecuilder/src/molecule_geometry.cpp
molecuilder/src/molecule_graph.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/tesselation.cpp
molecuilder/src/tesselationhelpers.cpp
molecuilder/src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp
molecuilder/src/unittests/bondgraphunittest.cpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/molecule_dynamics.cpp

    r90c4460 r0d111b  
    1414#include "molecule.hpp"
    1515#include "parser.hpp"
     16#include "Plane.hpp"
    1617
    1718/************************************* Functions for class molecule *********************************/
    … …  
    3839    // determine normalized trajectories direction vector (n1, n2)
    3940    Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    40     trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
    41     trajectory1.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
     41    trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);
    4242    trajectory1.Normalize();
    4343    Norm1 = trajectory1.Norm();
    4444    Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    45     trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
    46     trajectory2.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));
     45    trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
    4746    trajectory2.Normalize();
    4847    Norm2 = trajectory1.Norm();
    4948    // check whether either is zero()
    5049    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    51       tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
     50      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.startstep));
    5251    } else if (Norm1 < MYEPSILON) {
    5352      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
    54       trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy first offset
    55       trajectory1.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));  // subtract second offset
    56       trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    57       trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
     53      trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Runner->Trajectory.R.at(Params.startstep);
     54      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
     55      trajectory1 -= trajectory2;   // project the part in norm direction away
    5856      tmp = trajectory1.Norm();  // remaining norm is distance
    5957    } else if (Norm2 < MYEPSILON) {
    6058      Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
    61       trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy second offset
    62       trajectory2.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));  // subtract first offset
    63       trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    64       trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
     59      trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);  // copy second offset
     60      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
     61      trajectory2 -= trajectory1;   // project the part in norm direction away
    6562      tmp = trajectory2.Norm();  // remaining norm is distance
    66     } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
     63    } else if ((fabs(trajectory1.ScalarProduct(trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
    6764  //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
    6865  //        Log() << Verbose(0) << trajectory1;
    6966  //        Log() << Verbose(0) << " and ";
    7067  //        Log() << Verbose(0) << trajectory2;
    71       tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
     68      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.startstep));
    7269  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    7370    } else { // determine distance by finding minimum distance
    … …  
    7976  //        Log() << Verbose(0) << trajectory2 << endl;
    8077      // determine normal vector for both
    81       normal.MakeNormalVector(&trajectory1, &trajectory2);
     78      normal = Plane(trajectory1, trajectory2,0).getNormal();
    8279      // print all vectors for debugging
    8380  //        Log() << Verbose(0) << "Normal vector in between: ";
    … …  
    8582      // setup matrix
    8683      for (int i=NDIM;i--;) {
    87         gsl_matrix_set(A, 0, i, trajectory1.x[i]);
    88         gsl_matrix_set(A, 1, i, trajectory2.x[i]);
    89         gsl_matrix_set(A, 2, i, normal.x[i]);
    90         gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep).x[i] - Runner->Trajectory.R.at(Params.startstep).x[i]));
     84        gsl_matrix_set(A, 0, i, trajectory1[i]);
     85        gsl_matrix_set(A, 1, i, trajectory2[i]);
     86        gsl_matrix_set(A, 2, i, normal[i]);
     87        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - Runner->Trajectory.R.at(Params.startstep)[i]));
    9188      }
    9289      // solve the linear system by Householder transformations
    … …  
    9693  //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    9794      // test whether we really have the intersection (by checking on c_1 and c_2)
    98       TestVector.CopyVector(&Runner->Trajectory.R.at(Params.startstep));
     95      trajectory1.Scale(gsl_vector_get(x,0));
    9996      trajectory2.Scale(gsl_vector_get(x,1));
    100       TestVector.AddVector(&trajectory2);
    10197      normal.Scale(gsl_vector_get(x,2));
    102       TestVector.AddVector(&normal);
    103       TestVector.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
    104       trajectory1.Scale(gsl_vector_get(x,0));
    105       TestVector.SubtractVector(&trajectory1);
     98      TestVector = Runner->Trajectory.R.at(Params.startstep) + trajectory2 + normal
     99                   - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
    106100      if (TestVector.Norm() < MYEPSILON) {
    107101  //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
    … …  
    172166    // first term: distance to target
    173167    Runner = Params.PermutationMap[Walker->nr];   // find target point
    174     tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)));
     168    tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.endstep)));
    175169    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    176170    result += Params.PenaltyConstants[0] * tmp;
    … …  
    231225    while(Runner->next != mol->end) {
    232226      Runner = Runner->next;
    233       Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)), Runner) );
     227      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(Runner->Trajectory.R.at(Params.endstep)), Runner) );
    234228    }
    235229  }
    … …  
    514508      Sprinter = mol->AddCopyAtom(Walker);
    515509      for (int n=NDIM;n--;) {
    516         Sprinter->x.x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
     510        Sprinter->x[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
    517511        // add to Trajectories
    518512        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
    519513        if (step < MaxSteps) {
    520           Walker->Trajectory.R.at(step).x[n] = Walker->Trajectory.R.at(startstep).x[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep).x[n] - Walker->Trajectory.R.at(startstep).x[n])*((double)step/(double)MaxSteps);
    521           Walker->Trajectory.U.at(step).x[n] = 0.;
    522           Walker->Trajectory.F.at(step).x[n] = 0.;
     514          Walker->Trajectory.R.at(step)[n] = Walker->Trajectory.R.at(startstep)[n] + (PermutationMap[Walker->nr]->Trajectory.R.at(endstep)[n] - Walker->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
     515          Walker->Trajectory.U.at(step)[n] = 0.;
     516          Walker->Trajectory.F.at(step)[n] = 0.;
    523517        }
    524518      }
    … …  
    582576    for(int i=0;i<AtomCount;i++)
    583577      for(int d=0;d<NDIM;d++) {
    584         Velocity.x[d] += Force.Matrix[0][i][d+5];
     578        Velocity[d] += Force.Matrix[0][i][d+5];
    585579      }
    586580    for(int i=0;i<AtomCount;i++)
    587581      for(int d=0;d<NDIM;d++) {
    588         Force.Matrix[0][i][d+5] -= Velocity.x[d]/(double)AtomCount;
     582        Force.Matrix[0][i][d+5] -= Velocity[d]/(double)AtomCount;
    589583      }
    590584    // solve a constrained potential if we are meant to
Note: See TracChangeset for help on using the changeset viewer.