Changeset f714979 for src


Ignore:
Timestamp:
Dec 8, 2008, 2:16:34 PM (16 years ago)
Author:
Christian Neuen <neuen@…>
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:
caf5d6
Parents:
a8bcea6
Message:

Another update w.r.t. the Tesselation.
Some signs switched, but atom indices might be misused.

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/boundary.cpp

    ra8bcea6 rf714979  
    514514  double a,b,c;
    515515
    516   Find_non_convex_border(out, tecplot, *TesselStruct, mol);
    517   /*
     516  //Find_non_convex_border(out, tecplot, *TesselStruct, mol); // Is now called from command line.
     517
    518518  // 1. calculate center of gravity
    519519  *out << endl;
    … …  
    559559
    560560  *out << Verbose(2) << "I created " << TesselStruct->TrianglesOnBoundaryCount << " triangles with " << TesselStruct->LinesOnBoundaryCount << " lines and " << TesselStruct->PointsOnBoundaryCount << " points." << endl;
    561   */
     561
    562562
    563563
    … …  
    582582  *out << Verbose(0) << "RESULT: The summed volume is " << setprecision(10) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3." << endl;
    583583
    584   /*
     584
    585585  // 7. translate all points back from CoG
    586586  *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
    … …  
    591591    Walker->x.Translate(CenterOfGravity);
    592592  }
    593   */
     593
    594594
    595595
    … …  
    10391039 *  with all neighbours of the candidate.
    10401040 */
    1041 void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
    1042 {
    1043   /* d2 is normal vector on the triangle
    1044    * d1 is normal on the triangle line, from which we come, as well as on d2.
     1041void Find_next_suitable_point(atom* a, atom* b, atom* Candidate, int n, Vector *Chord, Vector *d1, Vector *OldNormal, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule* mol, bool problem)
     1042{
     1043  /* OldNormal is normal vector on the old triangle
     1044   * d1 is normal on the triangle line, from which we come, as well as on OldNormal.
    10451045   */
    10461046  Vector dif_a; //Vector from a to candidate
    … …  
    10661066  {
    10671067
    1068           if (Chord->Norm()/(2*sin(dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find of Ball will touch atom
     1068          if (Chord->Norm()/(2*sin(0.5*dif_a.Angle(&dif_b)))<RADIUS) //Using Formula for relation of chord length with inner angle to find if Ball will touch atom
    10691069          {
    10701070                  if (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1))>Storage[0]) //This will give absolute preference to those in "right-hand" quadrants
    … …  
    10721072                          Opt_Candidate = Candidate;
    10731073                          Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
    1074                           Storage[1]=AngleCheck.Angle(d2);
     1074                          Storage[1]=AngleCheck.Angle(OldNormal);
    10751075                  }
    10761076                  else
    10771077                  {
    1078                           if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]< AngleCheck.Angle(d2)) or \
    1079                                           (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]> AngleCheck.Angle(d2)))
     1078                          if ((dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]>0 &&  Storage[1]> AngleCheck.Angle(OldNormal)) or \
     1079                                          (dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1)) == Storage[0] && Storage[0]<0 &&  Storage[1]< AngleCheck.Angle(OldNormal)))
    10801080                                  //Depending on quadrant we prefer higher or lower atom with respect to Triangle normal first.
    10811081                          {
    10821082                                  Opt_Candidate = Candidate;
    10831083                                  Storage[0]=dif_a.ScalarProduct(d1)/fabs(dif_a.ScalarProduct(d1));
    1084                                   Storage[1]=AngleCheck.Angle(d2);
     1084                                  Storage[1]=AngleCheck.Angle(OldNormal);
    10851085                          }
    10861086                          else
    10871087                          {
    10881088                                  if (problem)
    1089                                           cout << "Loses to better candidate" << endl;
     1089                                          cout << "Looses to better candidate" << endl;
    10901090                          }
    10911091                  }
    … …  
    11091109                  Walker = mol->ListOfBondsPerAtom[Candidate->nr][i]->GetOtherAtom(Candidate);
    11101110
    1111                   Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, d2, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
     1111                  Find_next_suitable_point(a, b, Walker, n+1, Chord, d1, OldNormal, Opt_Candidate, Storage, RADIUS, mol, problem);  //call function again
    11121112      }
    11131113  }
    … …  
    11171117 * this function fins a triangle to a line, adjacent to an existing one.
    11181118 */
    1119 void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS)
     1119void Tesselation::Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N)
    11201120{
    11211121        cout << Verbose(1) << "Looking for next suitable triangle \n";
    … …  
    11471147  direction1.VectorProduct(&(T.NormalVector));
    11481148
    1149   if (direction1.ScalarProduct(&helper)>0)
     1149  if (direction1.ScalarProduct(&helper)<0)
    11501150    {
    11511151      direction1.Scale(-1);
    … …  
    11581158  Find_next_suitable_point(Line.endpoints[0]->node, Line.endpoints[1]->node, Line.endpoints[0]->node, 0, &Chord, &direction1, &(T.NormalVector), Opt_Candidate, Storage, RADIUS, mol,0);
    11591159
    1160   if (Opt_Candidate==NULL)
     1160  if (N>-1)
    11611161  {
    11621162          cout << Verbose(1) << "No new Atom found, triangle construction will crash" << endl;
    … …  
    11801180  cout << " Optimal candidate is " << *Opt_Candidate << endl;
    11811181  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
    1182   cout << " Optimal candidate is " << *Opt_Candidate << endl;
     1182
    11831183  BPS[1] = new class BoundaryPointSet(Line.endpoints[1]->node);
    11841184  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    … …  
    11901190
    11911191  cout << Verbose(1) << "Checking for already present lines ... " << endl;
    1192   for(int i=0;i<NDIM;i++) // sind Linien bereits vorhanden ???
    1193     {
    1194 
    1195 
    1196   if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
    1197         {
    1198                 LinesOnBoundaryCount++;
    1199         }
    1200     }
     1192  for(int i=0; i<NDIM; i++) // sind Linien bereits vorhanden ???
     1193  {
     1194          if ( (LinesOnBoundary.insert( LinePair(LinesOnBoundaryCount, BTS->lines[i]))).second);
     1195          {
     1196                  LinesOnBoundaryCount++;
     1197          }
     1198  }
    12011199  cout << Verbose(1) << "Constructing normal vector for this triangle ... " << endl;
    12021200
    … …  
    12191217        atom* Walker;
    12201218
    1221         AngleCheck.CopyVector(&(Candidate->x));
    1222         AngleCheck.SubtractVector(&(a->x));
    1223         if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
    1224     {
    1225           //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
    1226       Opt_Candidate=Candidate;
    1227       Storage[0]=AngleCheck.Angle(&Oben);
    1228       //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
    1229     };
     1219        if (a->nr !=Candidate->nr)
     1220        {
     1221                AngleCheck.CopyVector(&(Candidate->x));
     1222                AngleCheck.SubtractVector(&(a->x));
     1223                if (AngleCheck.Norm() < RADIUS and AngleCheck.Angle(&Oben) < Storage[0])
     1224                {
     1225                        //cout << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
     1226                        Opt_Candidate=Candidate;
     1227                        Storage[0]=AngleCheck.Angle(&Oben);
     1228                        //cout << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[1]);
     1229                }
     1230                else{
     1231                        if (AngleCheck.Norm() > RADIUS)
     1232                        {
     1233                                cout << Verbose(1) << "refused due to Radius" << AngleCheck.Norm() << endl;
     1234                        }
     1235                        else{
     1236                                cout << Verbose(1) << "Supposedly looses to better candidate" << Opt_Candidate->nr << endl;
     1237                        }
     1238                }
     1239        }
     1240
    12301241        if (n<5)
    12311242    {
    … …  
    12461257  int i=0;
    12471258  atom* Walker;
    1248   atom* Walker2;
    1249   atom* Walker3;
     1259  atom* FirstPoint;
     1260  atom* SecondPoint;
    12501261  int max_index[3];
    12511262  double max_coordinate[3];
    … …  
    12671278    {
    12681279        Walker = Walker->next;
    1269       for (i=0; i<3; i++)
    1270       {
    1271           if (Walker->x.x[i] > max_coordinate[i])
    1272             {
    1273               max_coordinate[i]=Walker->x.x[i];
    1274               max_index[i]=Walker->nr;
    1275             }
    1276       }
    1277     }
    1278 
    1279 cout << Verbose(1) << "Found starting atom \n";
     1280        for (i=0; i<3; i++)
     1281        {
     1282                if (Walker->x.x[i] > max_coordinate[i])
     1283                {
     1284                        max_coordinate[i]=Walker->x.x[i];
     1285                        max_index[i]=Walker->nr;
     1286                }
     1287        }
     1288    }
     1289
     1290  cout << Verbose(1) << "Found maximum coordinates. "<< endl;
    12801291  //Koennen dies fuer alle Richtungen, legen hier erstmal Richtung auf k=0
    1281   const int k=2;
     1292  const int k=0;
    12821293
    12831294  Oben.x[k]=1.;
    1284   Walker = mol->start;
    1285   Walker = Walker->next;
    1286   while (Walker->nr != max_index[k])
     1295  FirstPoint = mol->start;
     1296  FirstPoint = FirstPoint->next;
     1297  while (FirstPoint->nr != max_index[k])
    12871298    {
    1288       Walker = Walker->next;
    1289     }
    1290 cout << Verbose(1) << "Number of start atom: " << Walker->nr << endl;
     1299          FirstPoint = FirstPoint->next;
     1300    }
     1301  cout << Verbose(1) << "Coordinates of start atom: " << FirstPoint->x.x[0] << endl;
    12911302  double Storage[2];
    12921303  atom* Opt_Candidate = NULL;
    12931304  Storage[0]=999999.;   // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
    12941305  Storage[1]=999999.;  // This will be an angle looking for the third point.
    1295 cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[Walker->nr] << endl;
    1296 
    1297   for (i=1; i< (mol->NumberOfBondsPerAtom[Walker->nr]); i++)
    1298     {
    1299           Walker2 = mol->ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
    1300       Find_second_point_for_Tesselation(Walker, Walker2, 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
    1301     }
    1302 
    1303   Walker2 = Opt_Candidate;
     1306  cout << Verbose(1) << "Number of Bonds: " << mol->NumberOfBondsPerAtom[FirstPoint->nr] << endl;
     1307
     1308  Find_second_point_for_Tesselation(FirstPoint, mol->ListOfBondsPerAtom[FirstPoint->nr][0]->GetOtherAtom(FirstPoint), 0, Oben, Opt_Candidate, Storage, mol, RADIUS);
     1309
     1310
     1311  SecondPoint = Opt_Candidate;
    13041312  Opt_Candidate=NULL;
    1305   helper.CopyVector(&(Walker->x));
    1306   helper.SubtractVector(&(Walker2->x));
     1313  helper.CopyVector(&(FirstPoint->x));
     1314  helper.SubtractVector(&(SecondPoint->x));
    13071315  Oben.ProjectOntoPlane(&helper);
    13081316  helper.VectorProduct(&Oben);
    … …  
    13101318  Storage[1]= 9999999.; // This will be an angle looking for the third point.
    13111319
    1312   Chord.CopyVector(&(Walker->x)); // bring into calling function
    1313   Chord.SubtractVector(&(Walker2->x));
     1320  Chord.CopyVector(&(FirstPoint->x)); // bring into calling function
     1321  Chord.SubtractVector(&(SecondPoint->x));
    13141322
    13151323  cout << Verbose(1) << "Looking for third point candidates \n";
    1316   Find_next_suitable_point(Walker, Walker2, mol->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
    1317 
    1318 
    1319   //Starting Triangle is Walker, Walker2, Opt_Candidate
    1320 
    1321   AddPoint(Walker);
    1322   AddPoint(Walker2);
     1324  Find_next_suitable_point(FirstPoint, SecondPoint, mol->ListOfBondsPerAtom[FirstPoint->nr][0]->GetOtherAtom(FirstPoint), 0, &Chord, &helper, &Oben, Opt_Candidate, Storage, RADIUS, mol, 0);
     1325
     1326
     1327  //Starting Triangle is Walker, SecondPoint, Opt_Candidate
     1328
     1329  AddPoint(FirstPoint);
     1330  AddPoint(SecondPoint);
    13231331  AddPoint(Opt_Candidate);
    13241332
    1325   cout << Verbose(1) << "The found starting triangle consists of " << *Walker << ", " << *Walker2 << " and " << *Opt_Candidate << "." << endl;
    1326 
    1327   BPS[0] = new class BoundaryPointSet(Walker);
    1328   BPS[1] = new class BoundaryPointSet(Walker2);
     1333  cout << Verbose(1) << "The found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " << *Opt_Candidate << "." << endl;
     1334
     1335  BPS[0] = new class BoundaryPointSet(FirstPoint);
     1336  BPS[1] = new class BoundaryPointSet(SecondPoint);
    13291337  BLS[0] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    1330   BPS[0] = new class BoundaryPointSet(Walker);
     1338  BPS[0] = new class BoundaryPointSet(FirstPoint);
    13311339  BPS[1] = new class BoundaryPointSet(Opt_Candidate);
    13321340  BLS[1] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    1333   BPS[0] = new class BoundaryPointSet(Walker);
    1334   BPS[1] = new class BoundaryPointSet(Walker2);
     1341  BPS[0] = new class BoundaryPointSet(Opt_Candidate);
     1342  BPS[1] = new class BoundaryPointSet(SecondPoint);
    13351343  BLS[2] = new class BoundaryLineSet(BPS , LinesOnBoundaryCount);
    13361344
    … …  
    13541362
    13551363
    1356 void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol)
    1357 {
     1364void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol)
     1365{
     1366        int N =0;
     1367        struct Tesselation *Tess = new Tesselation;
    13581368        cout << Verbose(1) << "Entering finding of non convex hull. " << endl;
    13591369        cout << flush;
    1360   const double RADIUS =6;
    1361   LineMap::iterator baseline;
    1362   Tess.Find_starting_triangle(mol, RADIUS);
    1363 
    1364   baseline = Tess.LinesOnBoundary.begin();
    1365   while (baseline != Tess.LinesOnBoundary.end()) {
    1366     if (baseline->second->TrianglesCount == 1)
    1367       {
    1368         cout << Verbose(1) << "Begin of Tesselation ... " << endl;
    1369                 Tess.Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS); //the line is there, so there is a triangle, but only one.
    1370         cout << Verbose(1) << "End of Tesselation ... " << endl;
    1371       }
    1372     else
    1373       {
    1374         cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
    1375       }
    1376   baseline++;
    1377   }
    1378 
    1379 };
     1370        const double RADIUS =6.;
     1371        LineMap::iterator baseline;
     1372        Tess->Find_starting_triangle(mol, RADIUS);
     1373
     1374        baseline = Tess->LinesOnBoundary.begin();
     1375        while (baseline != Tess->LinesOnBoundary.end()) {
     1376                if (baseline->second->TrianglesCount == 1)
     1377                {
     1378                        cout << Verbose(1) << "Begin of Tesselation ... " << endl;
     1379                        Tess->Find_next_suitable_triangle(out, tecplot, mol, *(baseline->second), *(baseline->second->triangles.begin()->second), RADIUS, N); //the line is there, so there is a triangle, but only one.
     1380                        cout << Verbose(1) << "End of Tesselation ... " << endl;
     1381                }
     1382                else
     1383                {
     1384                        cout << Verbose(1) << "There is a line with " << baseline->second->TrianglesCount << " triangles adjacent";
     1385                }
     1386                N++;
     1387                baseline++;
     1388        }
     1389
     1390};

  • src/boundary.hpp

    ra8bcea6 rf714979  
    8383    void AddPoint(atom * Walker);
    8484    void Find_starting_triangle(molecule* mol, const double RADIUS);
    85     void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS);
     85    void Find_next_suitable_triangle(ofstream *out, ofstream *tecplot, molecule* mol, BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, int N);
    8686
    8787    PointMap PointsOnBoundary;
    … …  
    106106void PrepareClustersinWater(ofstream *out, config *configuration, molecule *mol, double ClusterVolume, double celldensity);
    107107void Find_next_suitable_point(atom a, atom b, atom Candidate, int n, Vector *d1, Vector *d2, atom*& Opt_Candidate, double *Storage, const double RADIUS, molecule *mol, bool problem);
    108 void Find_non_convex_border(ofstream *out, ofstream *tecplot, Tesselation Tess, molecule* mol);
     108void Find_non_convex_border(ofstream *out, ofstream *tecplot, molecule* mol);
    109109
    110110

  • src/builder.cpp

    ra8bcea6 rf714979  
    801801            cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
    802802            cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
    803             cout << "\t-A <source>\tCreate adjacenzy list from bonds parsed from 'dbond'-style file." <<endl;
     803            cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
    804804            cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
    805805            cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    … …  
    813813            cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    814814            cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
     815            cout << "\t-N\tGet non-convex-envelope." << endl;
    815816            cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
    816817            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
    … …  
    10151016                        input->close();
    10161017                }
     1018                break;
     1019            case 'N':
     1020                ExitFlag = 1;
     1021                if ((argptr >= argc) || (argv[argptr][0] == '-')){
     1022                        ExitFlag = 255;
     1023                        cerr << "Not enough or invalid arguments given for non-convex envelope: -o <tecplot output file>" << endl;
     1024                        }
     1025                else {
     1026                        cout << Verbose(0) << "Evaluating npn-convex envelope.";
     1027                        ofstream *output = new ofstream(argv[argptr], ios::trunc);
     1028                        cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
     1029                        Find_non_convex_border((ofstream *)&cout, output, mol);
     1030                        output->close();
     1031                        delete(output);
     1032                        argptr+=1;
     1033                        }
    10171034                break;
    10181035            case 'T':

  • src/molecules.cpp

    ra8bcea6 rf714979  
    16771677
    16781678/** Creates an adjacency list of the molecule.
     1679 * We obtain an outside file with the indices of atoms which are bondmembers.
     1680 */
     1681void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)
     1682{
     1683
     1684        // 1 We will parse bonds out of the dbond file created by tremolo.
     1685                        int atom1, atom2, temp;
     1686                        atom *Walker, *OtherWalker;
     1687
     1688                if (!input)
     1689                {
     1690                        cout << Verbose(1) << "Opening silica failed \n";
     1691                };
     1692
     1693                        *input >> ws >> atom1;
     1694                        *input >> ws >> atom2;
     1695                cout << Verbose(1) << "Scanning file\n";
     1696                while (!input->eof()) // Check whether we read everything already
     1697                {
     1698                                *input >> ws >> atom1;
     1699                                *input >> ws >> atom2;
     1700                        if(atom2<atom1) //Sort indices of atoms in order
     1701                        {
     1702                                temp=atom1;
     1703                                atom1=atom2;
     1704                                atom2=temp;
     1705                        };
     1706
     1707                        Walker=start;
     1708                        while(Walker-> nr != atom1) // Find atom corresponding to first index
     1709                        {
     1710                                Walker = Walker->next;
     1711                        };
     1712                        OtherWalker = Walker->next;
     1713                        while(OtherWalker->nr != atom2) // Find atom corresponding to second index
     1714                        {
     1715                                OtherWalker= OtherWalker->next;
     1716                        };
     1717                        AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
     1718                        BondCount++; //Increase bond count of molecule
     1719                }
     1720
     1721                CreateListOfBondsPerAtom(out);
     1722
     1723};
     1724
     1725
     1726/** Creates an adjacency list of the molecule.
    16791727 * Generally, we use the CSD approach to bond recognition, that is the the distance
    16801728 * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
    … …  
    16941742 * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
    16951743 */
    1696 void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)// ofstream* out, double bonddistance, bool IsAngstroem)
    1697 {
    1698 
    1699         // 1 We will parse bonds out of the dbond file created by tremolo.
    1700                         int atom1, atom2, temp;
    1701                         atom *Walker, *OtherWalker;
    1702 
    1703                 if (!input)
    1704                 {
    1705                         cout << Verbose(1) << "Opening silica failed \n";
    1706                 };
    1707 
    1708                         *input >> ws >> atom1;
    1709                         *input >> ws >> atom2;
    1710                 cout << Verbose(1) << "Scanning file\n";
    1711                 while (!input->eof()) // Check whether we read 2 items
    1712                 {
    1713                                 *input >> ws >> atom1;
    1714                                 *input >> ws >> atom2;
    1715                         if(atom2<atom1) //Sort indizes of atoms in order
    1716                         {
    1717                                 temp=atom1;
    1718                                 atom1=atom2;
    1719                                 atom2=temp;
    1720                         };
    1721 
    1722                         Walker=start;
    1723                         while(Walker-> nr != atom1) // Find atom corresponding to first index
    1724                         {
    1725                                 Walker = Walker->next;
    1726                         };
    1727                         OtherWalker = Walker->next;
    1728                         while(OtherWalker->nr != atom2) // Find atom corresponding to second index
    1729                         {
    1730                                 OtherWalker= OtherWalker->next;
    1731                         };
    1732                         AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
    1733                         BondCount++; //Increase bond count of molecule
    1734                 }
    1735 
    1736                 CreateListOfBondsPerAtom(out);
    1737 
    1738 };
    17391744void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
    17401745{
Note: See TracChangeset for help on using the changeset viewer.