Changeset cccac6


Ignore:
Timestamp:
Jun 22, 2008, 1:27:47 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
7a40c8
Parents:
b73d28 (diff), 73bc6b (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 'bisect' into AdaptiveMolecuilder

Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    rb73d28 rcccac6  
    11# exclude some dirs
    2 build64
     2build*
    33bin
    44tests
    … …  
    66# exclude generated html
    77*.html
     8
     9# exclude Eclipse files
     10.cdtproject
     11.cproject
     12.project
     13.settings
     14 
    815
    916# exclude automake generated stuff

  • molecuilder/src/parser.cpp

    rb73d28 rcccac6  
    186186      for(int k=skipcolumns;k--;)
    187187        lines >> filename;
    188       for(int k=0;(k<ColumnCounter) && (!line.eof());k++) {
     188      for(int k=0;(k<ColumnCounter) && (!lines.eof());k++) {
    189189        lines >> Matrix[i][j][k];
    190190        //cout << " " << setprecision(2) << Matrix[i][j][k];;

  • pcp/src/output.c

    rb73d28 rcccac6  
    316316{
    317317  int i,j,k, Index, owner;
    318   size_t zahl;
     318  int zahl;
    319319  struct Lattice *Lat = &P->Lat;
    320320  //struct RunStruct *R = &P->R;
    … …  
    577577  int LevelNo, readnr=0;
    578578  int breaksignal = test ? 1 : 2; // 0 - ok, 1 - test failed, 2 - throw Error
    579   size_t zahl;
     579  int zahl;
    580580  char suffixdat[MAXSTRINGSIZE], suffixdoc[MAXSTRINGSIZE];
    581581  int Num = 0, colorNo = 0;

  • util/src/NanoCreator.c

    rb73d28 rcccac6  
    55#include <time.h>
    66
     7#include <gsl/gsl_matrix.h>
     8#include <gsl/gsl_vector.h>
     9#include <gsl/gsl_eigen.h>
     10#include <gsl/gsl_blas.h>
     11
    712#include "NanoCreator.h"
     13
    814
    915// ---------------------------------- F U N C T I O N S ----------------------------------------------
    … …  
    404410    exit(255);
    405411  }
    406   double betrag1, betrag2;
     412  double betrag;
    407413  int i;
    408414 
    409415  // first vector is untouched
    410416  // second vector
    411   betrag1 = Projection(orthvector[axis[1]], orthvector[axis[0]]);
     417  betrag = Projection(orthvector[axis[1]], orthvector[axis[0]]);
     418  fprintf(stdout,"%lg\t",betrag);
    412419  for (i=0;i<NDIM;i++)
    413      orthvector[axis[0]][i] -= orthvector[axis[1]][i] * betrag1;
     420     orthvector[axis[1]][i] -= orthvector[axis[0]][i] * betrag;
    414421  // third vector
    415   betrag1 = Projection(orthvector[axis[0]], orthvector[axis[2]]);
    416   betrag2 = Projection(orthvector[axis[1]], orthvector[axis[2]]);
     422  betrag = Projection(orthvector[axis[0]], orthvector[axis[2]]);
     423  fprintf(stdout,"%lg\t",betrag);
    417424  for (i=0;i<NDIM;i++)
    418      orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag1 + orthvector[axis[1]][i] * betrag2;
     425     orthvector[axis[2]][i] -= orthvector[axis[0]][i] * betrag;
     426  betrag = Projection(orthvector[axis[1]], orthvector[axis[2]]);
     427  fprintf(stdout,"%lg\n",betrag);
     428  for (i=0;i<NDIM;i++)
     429     orthvector[axis[2]][i] -= orthvector[axis[1]][i] * betrag;
    419430}
    420431
    … …  
    581592}
    582593
     594/** Creates orthogonal vectors to directions axis[0] and axis[1], assuming that axis[2] is always orthogonal.
     595 * \param **OrthoVector contains vectors set on return with axis[2] equal to Vector[axis[2]]
     596 * \param **Vector vectors to orthogonalize against
     597 * \param *axis lookup for which direction is which.
     598 */
     599void CreateOrthogonalAxisVectors(double **OrthoVector, double **Vector, int *axis) {
     600  int i,j;
     601  double factor;
     602  // allocate memory
     603  int *TempAxis = (int *) Malloc(sizeof(int)*NDIM, "Main: *TempAxis");
     604  double **TempVectors = (double **) Malloc(sizeof(double *)*NDIM, "Main: *TempVectors");
     605  for (i=0; i<NDIM; i++)
     606    TempVectors[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: TempVectors");
     607 
     608  for (i=0; i<NDIM; i++)
     609    for (j=0; j<NDIM; j++)
     610      TempVectors[i][j] = Vector[i][j];
     611  // GramSchmidt generates Vector[1] orthogonal to Vector[0] and Vector[2] ortho. to Vector[1] and Vector[0]!
     612  TempAxis[0] = axis[2];  // (axis 2 is the orthogonal plane axis)
     613  TempAxis[1] = axis[0];
     614  TempAxis[2] = axis[1];
     615  Orthogonalize(TempVectors, TempAxis);
     616  factor = Norm(Vector[axis[0]])/Norm(TempVectors[TempAxis[2]]);
     617  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[0]]) > 0) ? 1. : -1.;
     618  for (i=0; i<NDIM; i++)
     619    OrthoVector[axis[0]][i] = TempVectors[TempAxis[2]][i]*factor;
     620 
     621  TempAxis[1] = axis[1];
     622  TempAxis[2] = axis[0];
     623  for (i=0; i<NDIM; i++)
     624    for (j=0; j<NDIM; j++)
     625      TempVectors[i][j] = Vector[i][j];
     626  Orthogonalize(TempVectors, TempAxis);
     627  factor = Norm(Vector[axis[1]])/Norm(TempVectors[TempAxis[2]]);
     628  factor *= (Projection(TempVectors[TempAxis[2]], Vector[axis[1]]) > 0) ? 1. : -1.;
     629  for (i=0; i<NDIM; i++)
     630    OrthoVector[axis[1]][i] = TempVectors[TempAxis[2]][i]*factor;
     631 
     632  for (i=0; i<NDIM; i++)
     633    OrthoVector[axis[2]][i] = Vector[axis[2]][i];
     634  //print vectors
     635  fprintf(stdout, "Orthogonal vectors are: \n");
     636  for (i=0; i<NDIM; i++) {
     637    for (j=0; j<NDIM; j++)
     638      fprintf(stdout, "%lg\t", OrthoVector[axis[i]][j]);
     639    fprintf(stdout, "\n");
     640  }
     641 
     642  // free memory
     643  Free(TempAxis, "CreateOrthogonalAxisVectors: *TempAxis");
     644  for (i=0; i<NDIM; i++ )
     645    Free(TempVectors[i], "CreateOrthogonalAxisVectors: TempVectors");
     646  Free(TempVectors, "CreateOrthogonalAxisVectors: *TempVectors");
     647};
     648
    583649// ================================= other functions ==============================
    584650
     651/** Prints a debug message.
     652 * \param *msg debug message
     653 */
    585654void Debug(char *msg)
    586655{
    … …  
    591660  fprintf(stdout, "%s", msg);
    592661}
     662
    593663
    594664// --------------------------------------- M A I N ---------------------------------------------------
    … …  
    597667  char *CellFilename = NULL, *SheetFilename = NULL, *TubeFilename = NULL, *TorusFilename = NULL;
    598668  char *SheetFilenameAligned = NULL, *TubeFilenameAligned = NULL;
    599         double **Vector, **Recivector = NULL, **Tubevector = NULL, **TubevectorInverse = NULL;
     669        double **Vector, **OrthoVector, **Recivector = NULL, **Tubevector = NULL, **TubevectorInverse = NULL;
    600670  double *atom = NULL, *atom_transformed = NULL;
    601671  double *x = NULL, *coord = NULL, *tempvector = NULL, *offset = NULL;
    … …  
    645715  atom = (double *) Malloc(sizeof(double)*NDIM, "Main: atom");
    646716  Vector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Vector");
     717  OrthoVector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *OrthoVector");
    647718  Recivector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Recivector");
    648719  Tubevector = (double **) Malloc(sizeof(double *)*NDIM, "Main: *Tubevector");
    … …  
    650721  for (i=0; i<NDIM; i++ )  {
    651722    Vector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Vector");
     723    OrthoVector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: OrthoVector");
    652724    Recivector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Recivector");
    653725    Tubevector[i] = (double *) Malloc(sizeof(double)*NDIM, "Main: Tubevector");
    … …  
    701773  PrintMatrix(stdout, Recivector);
    702774  fprintf(stdout, "Reciprocal volume is %lg\n", Determinant(Recivector));
    703 
     775 
    704776  fprintf(stdout, "Vector magnitudes: %5.5lg %5.5lg %5.5lg\n", Norm(Vector[0]), Norm(Vector[1]), Norm(Vector[2]));
    705777   
    … …  
    788860    fscanf(stdin, "%d %d %d", &axis[0], &axis[1], &axis[2]);
    789861    fprintf(stdout, "axis: %d %d %d\n", axis[0], axis[1], axis[2]);
     862   
     863    // create orthogonal vectors individually for each unit cell vector
     864    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
     865    fprintf(stdout, "Orthogonal vector axis[0] %lg\n", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
     866    fprintf(stdout, "Orthogonal vector axis[1] %lg\n", Projection(Vector[axis[1]], OrthoVector[axis[1]]));
     867
    790868    do {     
    791869      fprintf(stdout, "\nNow specify the two natural numbers (m n) defining the chiral angle, \nif the result is crap, try flipping to (m,n): ");
    … …  
    795873      fprintf(stdout, "chiral0: %d\tchiral1: %d\n", chiral[0], chiral[1]);
    796874      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];
     875        Tubevector[axis[0]][i] = (double)chiral[0] * Vector[axis[0]][i] + (double)chiral[1] * Vector[axis[1]][i];
     876        //Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
    799877        //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];
     878        //Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     879        Tubevector[axis[1]][i] = (double)chiral[0] * OrthoVector[axis[0]][i] - (double)chiral[1] * OrthoVector[axis[1]][i];
    801880        //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];
    802881//        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",
    … …  
    806885        Tubevector[axis[2]][i] = Vector[axis[2]][i];
    807886      }
     887      // here we assume, that Vector[axis[2]] is along z direction!
     888      gsl_matrix *M = gsl_matrix_alloc(2,2);
     889      gsl_matrix *C = gsl_matrix_alloc(2,2);
     890      gsl_matrix *evec = gsl_matrix_alloc(2,2);
     891      gsl_vector *eval = gsl_vector_alloc(2);
     892      gsl_vector *v = gsl_vector_alloc(2);
     893      gsl_vector *u = gsl_vector_alloc(2);
     894      gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
     895      gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
     896      gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
     897      gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
     898      gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
     899      gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
     900      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     901      gsl_eigen_symmv(M, eval, evec, w);
     902      gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
     903      fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
     904      fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
     905      gsl_matrix_set(M, 0,0, 0.);
     906      gsl_matrix_set(M, 1,0, 1.);
     907      gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
     908      gsl_matrix_set(M, 1,1, 0.);
     909      gsl_vector_set(v,0,(double)chiral[0]);
     910      gsl_vector_set(v,1,(double)chiral[1]);
     911      gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
     912      gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
     913      fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     914      gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     915      fprintf(stdout, "Looking for factor to integer...\n");
     916      for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
     917        x1 = gsl_vector_get(u,0)*(double)i;
     918        x2 = gsl_vector_get(u,1)*(double)i;
     919        x3 =
     920        fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
     921        if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-6) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-6 )) {
     922          gsl_blas_dscal((double)i, u);
     923          break;
     924        }
     925      }
     926      fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));
     927
     928      // get length
     929      double x[NDIM];
     930      for (i=0;i<NDIM;i++)
     931        x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     932      angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
     933      fprintf(stdout, "angle is %lg\n", angle);
     934      for (i=0;i<NDIM;i++) {
     935        Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     936      }
     937
     938      // Probe
     939      gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
     940      gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
     941      gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
     942      gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
     943      gsl_vector_set(v,0,(double)chiral[0]);
     944      gsl_vector_set(v,1,(double)chiral[1]);
     945      gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
     946      gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     947      x1=1.;
     948      gsl_blas_ddot(u,eval,&x1);
     949      fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);
     950
     951      gsl_matrix_free(M);
     952      gsl_matrix_free(C);
     953      gsl_matrix_free(evec);
     954      gsl_vector_free(eval);
     955      gsl_vector_free(v);
     956      gsl_vector_free(u);
     957      gsl_eigen_symmv_free(w);
     958     
     959      if (fabs(x1) > 1e-6) {
     960        fprintf(stderr,"Resulting TubeVectors of axis %d and %d and not orthogonal, aborting.\n", axis[0], axis[1]);
     961        return(128);
     962      }
     963       
     964     
     965      angle = Projection(Tubevector[axis[1]], Vector[axis[0]]);
     966      fprintf(stdout, "Projection Tubevector1 axis[0] %lg %lg\n", angle, 1./angle);
     967      angle = Projection(Tubevector[axis[1]], Vector[axis[1]]);
     968      fprintf(stdout, "Projection Tubevector1 axis[1] %lg %lg\n", angle, 1./angle);
     969     
    808970/*      fprintf(stdout, "Vector\n");
    809971      PrintMatrix(stdout, Vector);
    … …  
    8951057    fprintf(stdout, "%d %d %d\n", factors[0], factors[1], factors[2]);
    8961058
     1059    // create orthogonal vectors individually for each unit cell vector
     1060    CreateOrthogonalAxisVectors(OrthoVector, Vector, axis);
     1061    fprintf(stdout, "Orthogonal vector axis[0] %lg", Projection(Vector[axis[0]], OrthoVector[axis[0]]));
     1062    fprintf(stdout, "Orthogonal vector axis[1] %lg", Projection(Vector[axis[1]], OrthoVector[axis[1]]));
    8971063    // create Tubevectors
    8981064    for (i=0;i<NDIM;i++) {
    8991065      Tubevector[axis[0]][i] = chiral[0] * Vector[axis[0]][i] + chiral[1] * Vector[axis[1]][i];
    9001066      //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];
     1067      //Tubevector[axis[1]][i] = -chiral[1] * Vector[axis[0]][i] + chiral[0] * Vector[axis[1]][i];
     1068      Tubevector[axis[1]][i] = chiral[0] * OrthoVector[axis[0]][i] + chiral[1] * OrthoVector[axis[1]][i];
    9021069      //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];
    9031070      Tubevector[axis[2]][i] = Vector[axis[2]][i];
     1071    }
     1072    // here we assume, that Vector[axis[2]] is along z direction!
     1073    gsl_matrix *M = gsl_matrix_alloc(2,2);
     1074    gsl_matrix *C = gsl_matrix_alloc(2,2);
     1075    gsl_matrix *evec = gsl_matrix_alloc(2,2);
     1076    gsl_vector *eval = gsl_vector_alloc(2);
     1077    gsl_vector *v = gsl_vector_alloc(2);
     1078    gsl_vector *u = gsl_vector_alloc(2);
     1079    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(2);
     1080    gsl_matrix_set(C, 0,0, Vector[axis[0]][0]);
     1081    gsl_matrix_set(C, 1,0, Vector[axis[0]][1]);
     1082    gsl_matrix_set(C, 0,1, Vector[axis[1]][0]);
     1083    gsl_matrix_set(C, 1,1, Vector[axis[1]][1]);
     1084    gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0, C, C, 0.0, M);
     1085    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     1086    gsl_eigen_symmv(M, eval, evec, w);
     1087    gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_ABS_DESC);
     1088    fprintf(stdout, "Eigenvalues: %lg\t%lg\n", gsl_vector_get(eval,0), gsl_vector_get(eval,1));
     1089    fprintf(stdout, "Eigenvectors: \t%lg\t%lg\n\t\t%lg\t%lg\n", gsl_matrix_get(evec,0,0), gsl_matrix_get(evec,0,1), gsl_matrix_get(evec,1,0), gsl_matrix_get(evec,1,1));
     1090    gsl_matrix_set(M, 0,0, 0.);
     1091    gsl_matrix_set(M, 1,0, 1.);
     1092    gsl_matrix_set(M, 0,1, -gsl_vector_get(eval,1)/gsl_vector_get(eval,0));
     1093    gsl_matrix_set(M, 1,1, 0.);
     1094    gsl_vector_set(v,0,(double)chiral[0]);
     1095    gsl_vector_set(v,1,(double)chiral[1]);
     1096    gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0, evec, M, 0.0, C);
     1097    gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0, C, evec, 0.0, M);
     1098    fprintf(stdout, "M: \t%lg\t%lg\n\t%lg\t%lg\n", gsl_matrix_get(M,0,0), gsl_matrix_get(M,0,1), gsl_matrix_get(M,1,0), gsl_matrix_get(M,1,1));
     1099    gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     1100    fprintf(stdout, "Looking for factor to integer...\n");
     1101    for(i=1;i<(chiral[0]+chiral[1])*(chiral[0]+chiral[1]);i++) {
     1102      x1 = gsl_vector_get(u,0)*(double)i;
     1103      x2 = gsl_vector_get(u,1)*(double)i;
     1104      x3 =
     1105      fprintf(stdout, "%d: %d\t%d vs. %lg\t%lg\n",i, ((int)(x1+x1/fabs(x1)*.5)), ((int)(x2+x2/fabs(x2)*.5)), (x1), (x2));
     1106      if (( fabs( ((int)(x1+x1/fabs(x1)*.5)) - (x1) ) < 1e-6) && ( fabs( ((int)(x2+x2/fabs(x2)*.5)) - (x2) ) < 1e-6 )) {
     1107        gsl_blas_dscal((double)i, u);
     1108        break;
     1109      }
     1110    }
     1111    fprintf(stdout, "(c,d) = (%lg,%lg)\n",gsl_vector_get(u,0), gsl_vector_get(u,1));
     1112
     1113    // get length
     1114    double x[NDIM];
     1115    for (i=0;i<NDIM;i++)
     1116      x[i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     1117    angle = Norm(x)/Norm(Tubevector[axis[1]]) ;//ScalarProduct(x,Tubevector[axis[1]])/Norm(Tubevector[axis[1]]);
     1118    fprintf(stdout, "angle is %lg\n", angle);
     1119    for (i=0;i<NDIM;i++) {
     1120      Tubevector[axis[1]][i] = gsl_vector_get(u,0) * Vector[axis[0]][i] + gsl_vector_get(u,1) * Vector[axis[1]][i];
     1121    }
     1122
     1123    // Probe
     1124    gsl_matrix_set(M, 0,0, Vector[axis[0]][0]);
     1125    gsl_matrix_set(M, 1,0, Vector[axis[0]][1]);
     1126    gsl_matrix_set(M, 0,1, Vector[axis[1]][0]);
     1127    gsl_matrix_set(M, 1,1, Vector[axis[1]][1]);
     1128    gsl_vector_set(v,0,(double)chiral[0]);
     1129    gsl_vector_set(v,1,(double)chiral[1]);
     1130    gsl_blas_dgemv(CblasNoTrans, 1.0, M, u, 0.0, eval);
     1131    gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, u);
     1132    x1=1.;
     1133    gsl_blas_ddot(u,eval,&x1);
     1134    fprintf(stdout, "Testing (c,d): (a,b) M^t M (c,d)^t = 0 ? : %lg\n", x1);
     1135
     1136    gsl_matrix_free(M);
     1137    gsl_matrix_free(C);
     1138    gsl_matrix_free(evec);
     1139    gsl_vector_free(eval);
     1140    gsl_vector_free(v);
     1141    gsl_vector_free(u);
     1142    gsl_eigen_symmv_free(w);   
     1143
     1144    if (fabs(x1) > 1e-6) {
     1145      fprintf(stderr,"Resulting TubeVectors of axis %d and %d and not orthogonal, aborting.\n", axis[0], axis[1]);
     1146      return(128);
    9041147    }
    9051148
    … …  
    9491192  }
    9501193
    951   //Orthogonalize(Tubevector,axis);
     1194  //int OrthOrder[NDIM] = { axis[2], axis[0], axis[1] };
     1195  //Orthogonalize(Tubevector,OrthOrder);
    9521196  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.);
     1197  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]]));
     1198  if ( angle/M_PI*180. > 90 ) {
     1199    fprintf(stderr, "There seems to be something wrong with the unit cell! for nanotube the angle should be 60 degrees for example!\n");
     1200    return 1;
     1201  }
    9541202  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.);
     1203  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]]));
    9561204  //angle -= acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
    9571205  //angle = 30./180.*M_PI - acos(Projection(Vector[axis[0]], Tubevector[axis[0]]));
    958   angle = - acos(Projection(Tubevector[axis[0]], Vector[axis[0]]));
     1206  angle = -acos(Projection(Tubevector[axis[0]], Vector[axis[0]]));
    9591207  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.);
    962   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",
     1208  if (fabs(Tubevector[axis[0]][0]) > MYEPSILON)
     1209    angle = M_PI/2. - acos(Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
     1210  else
     1211    angle = 0.;
     1212  fprintf(stdout, "The absolute alignment rotation angle then is %lg %lg\n", angle/M_PI*180., Tubevector[axis[0]][0]/Norm(Tubevector[axis[0]]));
     1213  fprintf(stdout, "\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. final torus radius of %5.5f A.\n",
    9631214    acos(Projection(Vector[axis[0]], Tubevector[axis[0]]))/M_PI*180.,
    9641215    (double)factors[0]*Norm(Tubevector[axis[0]])/(2.*M_PI),
    … …  
    9681219  Orthogonalize(Tubevector, axis);   // with correct translational vector, not needed anymore (? what's been done here. Hence, re-inserted)
    9691220  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");
     1221  fprintf(stdout, "Tubevectors are \n");
    9711222  PrintMatrix(stdout, Tubevector);
    9721223  MatrixInversion(Tubevector, TubevectorInverse);
    … …  
    10581309            // project down on major and minor Tubevectors and check for length if out of sheet
    10591310            tempvector = MatrixTrafoInverse(coord, TubevectorInverse);
    1060             if (((tempvector[axis[0]] + MYEPSILON) >= 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
    1061                 ((tempvector[axis[1]] + MYEPSILON) >= 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
    1062                 ((tempvector[axis[2]] + MYEPSILON) >= 0) && ((factors[2] - tempvector[axis[2]]) > MYEPSILON)) { // check if within rotated cell          numbersheet++;
     1311            if (((tempvector[axis[0]] + MYEPSILON) > 0) && ((factors[0] - tempvector[axis[0]]) > MYEPSILON) &&
     1312                ((tempvector[axis[1]] + MYEPSILON) > 0) && ((factors[1] - tempvector[axis[1]]) > MYEPSILON) &&
     1313                ((tempvector[axis[2]] + MYEPSILON) > 0) && ((factors[2] - tempvector[axis[2]]) > MYEPSILON)) { // check if within rotated cell          numbersheet++;
    10631314              //if (nummer >= 2)  strcpy(atombuffer[nr].name, "O");
    10641315              //nummer++;
    … …  
    11661417      fprintf(TubeFile, "%s\t%lg\t%lg\t%lg\n", name, atom_transformed[0], atom_transformed[1], atom_transformed[2]);
    11671418      // 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);
     1419      x1 = atom_transformed[0]*cos(-angle) + atom_transformed[1] * sin(-angle);
     1420      x2 = atom_transformed[0]*(-sin(-angle)) + atom_transformed[1] * cos(-angle);
    11701421      x3 = atom_transformed[2];
    1171       fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x1, x3, x2);
     1422      fprintf(TubeFileAligned, "%s\t%lg\t%lg\t%lg\n", name, x3, x2, x1);  // order so that symmetry is along z axis
    11721423      //fprintf(stdout, "%s\t%5.5lg\t%5.5lg\t%5.5lg\n", name, atom_transformed[0], atom_transformed[1] ,atom_transformed[2]);
    11731424
Note: See TracChangeset for help on using the changeset viewer.