Changeset 1b192d


Ignore:
Timestamp:
Jun 20, 2008, 5:55:00 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
e2f3fde
Parents:
c6a100
Message:

Newest tries to get NanoCreator finally working for all chiral angles. Nanotubes work for one direction, the other not. Hence, Nanotorus seam is clearly visible

File:
1 edited

Legend:

Unmodified
Added
Removed

  • util/src/NanoCreator.c

    rc6a100 r1b192d  
    404404    exit(255);
    405405  }
    406   double betrag1, betrag2;
     406  double betrag;
    407407  int i;
    408408 
    409409  // first vector is untouched
    410410  // second vector
    411   betrag1 = Projection(orthvector[axis[1]], orthvector[axis[0]]);
     411  betrag = Projection(orthvector[axis[1]], orthvector[axis[0]]);
     412  fprintf(stdout,"%lg\t",betrag);
    412413  for (i=0;i<NDIM;i++)
    413      orthvector[axis[0]][i] -= orthvector[axis[1]][i] * betrag1;
     414     orthvector[axis[1]][i] -= orthvector[axis[0]][i] * betrag;
    414415  // third vector
    415   betrag1 = Projection(orthvector[axis[0]], orthvector[axis[2]]);
    416   betrag2 = Projection(orthvector[axis[1]], orthvector[axis[2]]);
     416  betrag = Projection(orthvector[axis[0]], orthvector[axis[2]]);
     417  fprintf(stdout,"%lg\t",betrag);
    417418  for (i=0;i<NDIM;i++)
    418      orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag1 + orthvector[axis[1]][i] * betrag2;
     419     orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag;
     420  betrag = Projection(orthvector[axis[1]], orthvector[axis[2]]);
     421  fprintf(stdout,"%lg\n",betrag);
     422  for (i=0;i<NDIM;i++)
     423     orthvector[axis[2]][i] -= orthvector[axis[1]][i] * betrag;
    419424}
    420425
    … …  
    581586}
    582587
     588/** Creates orthogonal vectors to directions axis[0] and axis[1], assuming that axis[2] is always orthogonal.
     589 * \param **OrthoVector contains vectors set on return with axis[2] equal to Vector[axis[2]]
     590 * \param **Vector vectors to orthogonalize against
     591 * \param *axis lookup for which direction is which.
     592 */
     593void CreateOrthogonalAxisVectors(double **OrthoVector, double **Vector, int *axis) {
     594  int i,j;
     595  double factor;
     596  // allocate memory
     597  int *TempAxis = (int *) Malloc(sizeof(int)*NDIM, "Main: *TempAxis");
     598  double **TempVectors = (double **) Malloc(sizeof(double *)*NDIM, "Main: *TempVectors");
     599  for (i=0; i<NDIM; i++)
     600    TempVectors[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: TempVectors");
     601 
     602  for (i=0; i<NDIM; i++)
     603    for (j=0; j<NDIM; j++)
     604      TempVectors[i][j] = Vector[i][j];
     605  // GramSchmidt generates Vector[1] orthogonal to Vector[0] and Vector[2] ortho. to Vector[1] and Vector[0]!
     606  TempAxis[0] = axis[2];  // (axis 2 is the orthogonal plane axis)
     607  TempAxis[1] = axis[0];
     608  TempAxis[2] = axis[1];
     609  Orthogonalize(TempVectors, TempAxis);
     610  factor = Norm(Vector[axis[0]])/Norm(TempVectors[TempAxis[2]]);
     611  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[0]]) > 0) ? 1. : -1.;
     612  for (i=0; i<NDIM; i++)
     613    OrthoVector[axis[0]][i] = TempVectors[TempAxis[2]][i]*factor;
     614 
     615  TempAxis[1] = axis[1];
     616  TempAxis[2] = axis[0];
     617  for (i=0; i<NDIM; i++)
     618    for (j=0; j<NDIM; j++)
     619      TempVectors[i][j] = Vector[i][j];
     620  Orthogonalize(TempVectors, TempAxis);
     621  factor = Norm(Vector[axis[1]])/Norm(TempVectors[TempAxis[2]]);
     622  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[1]]) > 0) ? 1. : -1.;
     623  for (i=0; i<NDIM; i++)
     624    OrthoVector[axis[1]][i] = TempVectors[TempAxis[2]][i]*factor;
     625 
     626  for (i=0; i<NDIM; i++)
     627    OrthoVector[axis[2]][i] = Vector[axis[2]][i];
     628  //print vectors
     629  fprintf(stdout, "Orthogonal vectors are: \n");
     630  for (i=0; i<NDIM; i++) {
     631    for (j=0; j<NDIM; j++)
     632      fprintf(stdout, "%lg\t", OrthoVector[axis[i]][j]);
     633    fprintf(stdout, "\n");
     634  }
     635 
     636  // free memory
     637  Free(TempAxis, "CreateOrthogonalAxisVectors: *TempAxis");
     638  for (i=0; i<NDIM; i++ )
     639    Free(TempVectors[i], "CreateOrthogonalAxisVectors: TempVectors");
     640  Free(TempVectors, "CreateOrthogonalAxisVectors: *TempVectors");
     641};
     642
    583643// ================================= other functions ==============================
    584644
     645/** Prints a debug message.
     646 * \param *msg debug message
     647 */
    585648void Debug(char *msg)
    586649{
    … …  
    591654  fprintf(stdout, "%s", msg);
    592655}
     656
    593657
    594658// --------------------------------------- M A I N ---------------------------------------------------
    … …  
    597661  char *CellFilename = NULL, *SheetFilename = NULL, *TubeFilename = NULL, *TorusFilename = NULL;
    598662  char *SheetFilenameAligned = NULL, *TubeFilenameAligned = NULL;
    599         double **Vector, **Recivector = NULL, **Tubevector = NULL, **TubevectorInverse = NULL;
     663        double **Vector, **OrthoVector, **Recivector = NULL, **Tubevector = NULL, **TubevectorInverse = NULL;
    600664  double *atom = NULL, *atom_transformed = NULL;
    601665  double *x = NULL, *coord = NULL, *tempvector = NULL, *offset = NULL;
    … …  
    645709  atom = (double *) Malloc(sizeof(double)*NDIM, "Main: atom");
    646710  Vector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Vector");
     711  OrthoVector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *OrthoVector");
    647712  Recivector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Recivector");
    648713  Tubevector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Tubevector");
    … …  
    650715  for (i=0; i<NDIM; i++ )  {
    651716    Vector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Vector");
     717    OrthoVector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: OrthoVector");
    652718    Recivector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Recivector");
    653719    Tubevector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Tubevector");
    … …  
    701767  PrintMatrix(stdout, Recivector);
    702768  fprintf(stdout, "Reciprocal volume is %lg\n", Determinant(Recivector));
    703 
     769 
    704770  fprintf(stdout, "Vector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Vector[0]), Norm(Vector[1]), Norm(Vector[2]));
    705771   
    … …  
    788854    fscanf(stdin, "%d %d %d", &axis[0], &axis[1], &axis[2]);
    789855    fprintf(stdout, "axis: %d %d %d\n", axis[0], axis[1], axis[2]);
     856   
     857    // create orthogonal vectors individually for each unit cell vector
     858    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
     859    fprintf(stdout, "Orthogonal vector axis[0] %lg\n", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
     860    fprintf(stdout, "Orthogonal vector axis[1] %lg\n", Projection(Vector[axis[1]], OrthoVector[axis[1]]));
     861
    790862    do {     
    791863      fprintf(stdout, "\nNow specify the two natural numbers (m n) defining the chiral angle, \nif the result is crap, try flipping to (m,n): ");
    … …  
    795867      fprintf(stdout, "chiral0: %d\tchiral1: %d\n", chiral[0], chiral[1]);
    796868      for (i=0;i<NDIM;i++) {
    797         //Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
    798         Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
     869        Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
     870        //Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
    799871        //Tubevector[axis[0]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
    800         Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     872        //Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     873        Tubevector[axis[1]][i] = (double)chiral[0] * OrthoVector[axis[0]][i] - (double)chiral[1] * OrthoVector[axis[1]][i];
    801874        //Tubevector[axis[1]][i] = (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[0]][i] + (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[1]][i];
    802875//        fprintf(stderr, "Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i]\n = %lg * %lg + %lg * %lg = %lg + %lg = %lg\n\n",
    … …  
    806879        Tubevector[axis[2]][i] = Vector[axis[2]][i];
    807880      }
     881//      tmp =  Norm(Tubevector[axis[1]]);
     882//      for (i=0;i<NDIM;i++)
     883//        Tubevector[axis[1]][i] *= Norm(Vector[axis[1]])*sqrt(3)*sqrt(chiral[0]*chiral[0]+chiral[1]*chiral[1]+chiral[0]*chiral[1])/(double)ggT/tmp;
     884      fprintf(stdout, "Orthogonal Tubevector axis[0->1] %lg\n", Projection(Tubevector[axis[0]], Tubevector[axis[1]]));
     885      fprintf(stdout, "Orthogonal Tubevector axis[0->2] %lg\n", Projection(Tubevector[axis[0]], Tubevector[axis[2]]));
     886      fprintf(stdout, "Orthogonal Tubevector axis[1->2] %lg\n", Projection(Tubevector[axis[1]], Tubevector[axis[2]]));
     887     
    808888/*      fprintf(stdout, "Vector\n");
    809889      PrintMatrix(stdout, Vector);
    … …  
    895975    fprintf(stdout, "%d %d %d\n", factors[0], factors[1], factors[2]);
    896976
     977    // create orthogonal vectors individually for each unit cell vector
     978    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
     979    fprintf(stdout, "Orthogonal vector axis[0] %lg", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
     980    fprintf(stdout, "Orthogonal vector axis[1] %lg", Projection(Vector[axis[1]], OrthoVector[axis[1]]));
    897981    // create Tubevectors
    898982    for (i=0;i<NDIM;i++) {
    899983      Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
    900984      //Tubevector[axis[0]][i] = (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[0]][i] + (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[1]][i];
    901       Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     985      //Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     986      Tubevector[axis[1]][i] = chiral[0] * OrthoVector[axis[0]][i] + chiral[1] * OrthoVector[axis[1]][i];
    902987      //Tubevector[axis[1]][i] = (-chiral[0]-2.*chiral[1])/(double)ggT * Vector[axis[0]][i] + (2.*chiral[0]+chiral[1])/(double)ggT * Vector[axis[1]][i];
    903988      Tubevector[axis[2]][i] = Vector[axis[2]][i];
    904989    }
     990   
    905991
    906992    // retrieve seed ...
    … …  
    9491035  }
    9501036
    951   //Orthogonalize(Tubevector,axis);
     1037  //int OrthOrder[NDIM] = { axis[2], axis[0], axis[1] };
     1038  //Orthogonalize(Tubevector,OrthOrder);
    9521039  angle = acos(Projection(Vector[axis[0]], Vector[axis[1]])); // calcs angle between shanks in unit cell
    953   fprintf(stdout, "The basic angle between the two shanks of the unit cell is %lg\n", angle/M_PI*180.);
     1040  fprintf(stdout, "The basic angle between the two shanks of the unit cell is %lg %lg\n", angle/M_PI*180., Projection(Vector[axis[0]], Vector[axis[1]]));
     1041  if ( angle/M_PI*180. > 90 ) {
     1042    fprintf(stderr, "There seems to be something wrong with the unit cell! for nanotube the angle should be 60 degrees for example!\n");
     1043    return 1;
     1044  }
    9541045  angle = acos(Projection(Tubevector[axis[0]], Tubevector[axis[1]])); // calcs angle between shanks in unit cell
    955   fprintf(stdout, "The basic angle between the two shanks of the tube unit cell is %lg\n", angle/M_PI*180.);
     1046  fprintf(stdout, "The basic angle between the two shanks of the tube unit cell is %lg %lg\n", angle/M_PI*180., Projection(Tubevector[axis[0]], Tubevector[axis[1]]));
    9561047  //angle -= acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
    9571048  //angle = 30./180.*M_PI - acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
    958   angle = - acos(Projection(Tubevector[axis[0]], Vector[axis[0]]));
     1049  angle = -acos(Projection(Tubevector[axis[0]], Vector[axis[0]]));
    9591050  fprintf(stdout, "The relative alignment rotation angle then is %lg\n", angle/M_PI*180.);
    960   angle += acos(Projection(Vector[axis[0]], Vector[axis[1]]))/2.;
    961   fprintf(stdout, "The absolute alignment rotation angle then is %lg\n", angle/M_PI*180.);
     1051  if (fabs(Tubevector[axis[0]][0]) > MYEPSILON)
     1052    angle = M_PI/2. - acos(Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
     1053  else
     1054    angle = 0.;
     1055  fprintf(stdout, "The absolute alignment rotation angle then is %lg %lg\n", angle/M_PI*180., Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
    9621056  fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n",
    9631057    acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
    … …  
    9681062  Orthogonalize(Tubevector, axis);   // with correct translational vector, not needed anymore (? what's been done here. Hence, re-inserted)
    9691063  fprintf(stdout, "Tubevector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Tubevector[0]), Norm(Tubevector[1]), Norm(Tubevector[2]));
    970   fprintf(stdout, "Tubebvectors are \n");
     1064  fprintf(stdout, "Tubevectors are \n");
    9711065  PrintMatrix(stdout, Tubevector);
    9721066  MatrixInversion(Tubevector, TubevectorInverse);
    … …  
    10571151            //PrintVector(stdout, coord);
    10581152            // project down on major and minor Tubevectors and check for length if out of sheet
    1059             tempvector = MatrixTrafoInverse(coord, TubevectorInverse);
     1153            tempvector = MatrixTrafoInverse(offset, TubevectorInverse);
    10601154            if (((tempvector[axis[0]] + MYEPSILON) >= 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
    10611155                ((tempvector[axis[1]] + MYEPSILON) >= 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
    … …  
    11661260      fprintf(TubeFile, "%s\t%lg\t%lg\t%lg\n", name, atom_transformed[0], atom_transformed[1], atom_transformed[2]);
    11671261      // rotate and flip to align tube in z-direction
    1168       x1 = atom_transformed[0]*cos(angle) + atom_transformed[1] * sin(angle);
    1169       x2 = atom_transformed[0]*(-sin(angle)) + atom_transformed[1] * cos(angle);
     1262      x1 = atom_transformed[0]*cos(-angle) + atom_transformed[1] * sin(-angle);
     1263      x2 = atom_transformed[0]*(-sin(-angle)) + atom_transformed[1] * cos(-angle);
    11701264      x3 = atom_transformed[2];
    1171       fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x1, x3, x2);
     1265      fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x3, x2, x1);  // order so that symmetry is along z axis
    11721266      //fprintf(stdout, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", name, atom_transformed[0], atom_transformed[1] ,atom_transformed[2]);
    11731267
Note: See TracChangeset for help on using the changeset viewer.