Changeset f5586e for pcp


Ignore:
Timestamp:
Apr 22, 2008, 8:43:24 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
90c027
Parents:
cc46b0
Message:

RealBasisQ[] removed

RealBasisQ[] is again replaced by sqrt(RealBasisSQ[]) as we are going to switch to a pure matrix transformation for calculating whether points are still within the periodic cell and so forth instead of "hoping" the simulation box is rectangular.

Location:
pcp/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • pcp/src/data.h

    rcc46b0 rf5586e  
    563563  double RealBasis[NDIM_NDIM];  //!< Coefficients of the basis vectors
    564564  double RealBasisSQ[NDIM];     //!< squared Norm of each basis vector
    565   double RealBasisQ[NDIM];     //!< Norm of each basis vector
     565  //double RealBasisQ[NDIM];     //!< Norm of each basis vector
    566566  double InvBasis[NDIM_NDIM];   //!< Matrix-wise inverted basis vectors
    567567  double ReciBasis[NDIM_NDIM];  //!< Coefficients of the transposed(!), inverse basis "matrix" (i.e. reciprocal basis)

  • pcp/src/init.c

    rcc46b0 rf5586e  
    148148  for (i=0; i < NDIM_NDIM; i++) Lat->ReciBasis[i] *= factor;            // SM(Lat->ReciBasis, factor, NDIM_NDIM);
    149149  for (i=0; i < NDIM; i++) {
     150    Lat->RealBasisSQ[i] = RNORMSQ3(&(Lat->RealBasis[i*NDIM]));
    150151    Lat->ReciBasisSQ[i] = RNORMSQ3(&(Lat->ReciBasis[i*NDIM]));
    151     Lat->RealBasisSQ[i] = RNORMSQ3(&(Lat->RealBasis[i*NDIM]));
    152     Lat->RealBasisQ[i] = sqrt(Lat->RealBasisSQ[i]);
    153152  }
    154153        // Prepares LatticeLevels
    … …  
    279278      if (i%NDIM == NDIM-1) fprintf(stderr,"\n");
    280279    }
    281     fprintf(stderr,"(%i) Lat->RealbasisQ =  ",P->Par.me);
    282     for (i=0; i < NDIM; i++)
    283       fprintf(stderr,"%lg ",Lat->RealBasisQ[i]);
    284     fprintf(stderr,"\n");
    285280  }
    286281 

  • pcp/src/perturbed.c

    rcc46b0 rf5586e  
    15361536        i0 = n[2]*NUp[2]+N[2]*NUp[2]*(n[1]*NUp[1]+N[1]*NUp[1]*n0*NUp[0]);
    15371537
    1538 //        if (fabs(truedist(Lat,x[index[1]],Wcentre[index[1]],index[1])) >= borderstart * Lat->RealBasisQ[index[1]]/2.)
     1538//        if (fabs(truedist(Lat,x[index[1]],Wcentre[index[1]],index[1])) >= borderstart * sqrt(Lat->RealBasisSQ[index[1]])/2.)
    15391539//          if (max[index[1]] < sawtooth(Lat,truedist(Lat,x[index[1]],Wcentre[index[1]],index[1]),index[1]) * TempPsiR [i0]) {
    15401540//            max[index[1]] = sawtooth(Lat,truedist(Lat,x[index[1]],Wcentre[index[1]],index[1]),index[1]) * TempPsiR [i0]; 
    … …  
    15431543//          }
    15441544//
    1545 //        if (fabs(truedist(Lat,x[index[3]],Wcentre[index[3]],index[3])) >= borderstart * Lat->RealBasisQ[index[3]]/2.)
     1545//        if (fabs(truedist(Lat,x[index[3]],Wcentre[index[3]],index[3])) >= borderstart * sqrt(Lat->RealBasisSQ[index[3]])/2.)
    15461546//          if (max[index[3]] < sawtooth(Lat,truedist(Lat,x[index[3]],Wcentre[index[3]],index[3]),index[3]) * TempPsiR [i0]) {
    15471547//            max[index[3]] = sawtooth(Lat,truedist(Lat,x[index[3]],Wcentre[index[3]],index[3]),index[3]) * TempPsiR [i0]; 
    … …  
    18811881
    18821882/** Returns sawtooth shaped profile for position operator within cell.
    1883  * This is a mapping from -L/2...L/2 (L = Lattice#RealBasisQ) to -L/2 to L/2 with a smooth transition:
     1883 * This is a mapping from -L/2...L/2 (L = length of unit cell derived from Lattice#RealBasisSQ) to -L/2 to L/2 with a smooth transition:
    18841884 * \f[
    18851885 *  f(x): x \rightarrow \left \{
    … …  
    18901890 *    \end{array} \right \}
    18911891 * \f]
    1892  * \param *Lat pointer to Lattice structure for Lattice#RealBasisQ
     1892 * \param *Lat pointer to Lattice structure for Lattice#RealBasisSQ
    18931893 * \param L parameter x
    1894  * \param index component index for Lattice#RealBasisQ
     1894 * \param index component index for Lattice#RealBasisSQ
    18951895 */
    18961896inline double sawtooth(struct Lattice *Lat, double L, const int index)
    18971897{
    1898   double axis = Lat->RealBasisQ[index];
     1898  double axis = sqrt(Lat->RealBasisSQ[index]);
    18991899  double sawstart = Lat->SawtoothStart;
    19001900  double sawend = 1. - sawstart;
    … …  
    19591959 
    19601960  for (n=0;n<N[index];n++) {
    1961     x = (double)n/(double)N[index] * Lat->RealBasisQ[index];
    1962     //fprintf(stderr,"(%i) x %e\t Axis/2 %e\n",P->Par.me, x, Lat->RealBasisQ[index]/2. );     
    1963     x = MinImageConv(Lat, x, Lat->RealBasisQ[index]/2., index);
     1961    x = (double)n/(double)N[index] * sqrt(Lat->RealBasisSQ[index]);
     1962    //fprintf(stderr,"(%i) x %e\t Axis/2 %e\n",P->Par.me, x, sqrt(Lat->RealBasisSQ[index])/2. );     
     1963    x = MinImageConv(Lat, x, sqrt(Lat->RealBasisSQ[index])/2., index);
    19641964    fprintf(stderr,"%e\t%e\n", x, sawtooth(Lat,x,index)); 
    19651965  } 
    … …  
    19761976inline double MinImageConv(struct Lattice *Lat, const double R, const double r, const int index)
    19771977{
    1978   double axis = Lat->RealBasisQ[index];
     1978  double axis = sqrt(Lat->RealBasisSQ[index]);
    19791979  double result = 0.;
    19801980
    … …  
    30403040            RMat33Vec3(x, Lat->RealBasis, fac);
    30413041            i0 = n[2]+N[2]*(n[1]+N[1]*(n0));  // the index of current density must match LocalSizeR!
    3042             chi[in+dex*NDIM] += MinImageConv(Lat,x[cross_lookup[0]], Lat->RealBasisQ[cross_lookup[0]]/2.,cross_lookup[0]) * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]; // x[cross(in,0)], Lat->RealBasisQ[cross_lookup[0]]/2.
    3043             chi[in+dex*NDIM] -= MinImageConv(Lat,x[cross_lookup[2]], Lat->RealBasisQ[cross_lookup[2]]/2.,cross_lookup[2]) * CurrentDensity[dex*NDIM+cross_lookup[3]][i0]; // x[cross(in,2)], Lat->RealBasisQ[cross_lookup[2]]/2.
     3042            chi[in+dex*NDIM] += MinImageConv(Lat,x[cross_lookup[0]], sqrt(Lat->RealBasisSQ[cross_lookup[0]])/2.,cross_lookup[0]) * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]; // x[cross(in,0)], sqrt(Lat->RealBasisSQ[cross_lookup[0]])/2.
     3043            chi[in+dex*NDIM] -= MinImageConv(Lat,x[cross_lookup[2]], sqrt(Lat->RealBasisSQ[cross_lookup[2]])/2.,cross_lookup[2]) * CurrentDensity[dex*NDIM+cross_lookup[3]][i0]; // x[cross(in,2)], sqrt(Lat->RealBasisSQ[cross_lookup[2]])/2.
    30443044//            if (in == dex) field[in][i0] =
    3045 //                truedist(Lat,x[cross_lookup[0]], Lat->RealBasisQ[c[0]]/2.,cross_lookup[0]) * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]
    3046 //              - truedist(Lat,x[cross_lookup[2]], Lat->RealBasisQ[c[2]]/2.,cross_lookup[2]) * CurrentDensity[dex*NDIM+cross_lookup[3]][i0];
    3047             //fprintf(stderr,"(%i) temporary susceptiblity \\chi[%i][%i] += %e * %e = r[%i] * CurrDens[%i][%i] = %e\n",P->Par.me,in,dex,(double)n[cross_lookup[0]]/(double)N[cross_lookup[0]]*(Lat->RealBasisQ[cross_lookup[0]]),CurrentDensity[dex*NDIM+cross_lookup[1]][i0],cross_lookup[0],dex*NDIM+cross_lookup[1],i0,chi[in*NDIM+dex]);
     3045//                truedist(Lat,x[cross_lookup[0]], sqrt(Lat->RealBasisSQ[c[0]])/2.,cross_lookup[0]) * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]
     3046//              - truedist(Lat,x[cross_lookup[2]], sqrt(Lat->RealBasisSQ[c[2]])/2.,cross_lookup[2]) * CurrentDensity[dex*NDIM+cross_lookup[3]][i0];
     3047            //fprintf(stderr,"(%i) temporary susceptiblity \\chi[%i][%i] += %e * %e = r[%i] * CurrDens[%i][%i] = %e\n",P->Par.me,in,dex,(double)n[cross_lookup[0]]/(double)N[cross_lookup[0]]*(sqrt(Lat->RealBasisSQ[cross_lookup[0]])),CurrentDensity[dex*NDIM+cross_lookup[1]][i0],cross_lookup[0],dex*NDIM+cross_lookup[1],i0,chi[in*NDIM+dex]);
    30483048          }
    30493049      chi[in+dex*NDIM] *= mu0*discrete_factor/(2.*Lat->Volume); // integral factor
    … …  
    32353235          // read transformed sigma at core position and MPI_Allreduce
    32363236          for (g=0;g<NDIM;g++) {
    3237             n[0][g] = floor(I->I[it].R[NDIM*ion+g]/Lat->RealBasisQ[g]*(double)N[g]);
    3238             n[1][g] = ceil(I->I[it].R[NDIM*ion+g]/Lat->RealBasisQ[g]*(double)N[g]);
    3239             x[1][g] = I->I[it].R[NDIM*ion+g]/Lat->RealBasisQ[g]*(double)N[g] - (double)n[0][g];
     3237            n[0][g] = floor(I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g]);
     3238            n[1][g] = ceil(I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g]);
     3239            x[1][g] = I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g] - (double)n[0][g];
    32403240            x[0][g] = 1. - x[1][g];
    32413241            //fprintf(stderr,"(%i) n_floor[%i] = %i\tn_ceil[%i] = %i --- x_floor[%i] = %e\tx_ceil[%i] = %e\n",P->Par.me, g,n[0][g], g,n[1][g], g,x[0][g], g,x[1][g]);
    … …  
    37443744        k = cross(i,1);
    37453745        if (GA[g].G[k] == GA[g_bar].G[k]) {
    3746           //b = truedist(Lat, Lat->RealBasisQ[j], Wcentre[j], j);
    3747           b = Lat->RealBasisQ[j];
     3746          //b = truedist(Lat, sqrt(Lat->RealBasisSQ[j]), Wcentre[j], j);
     3747          b = sqrt(Lat->RealBasisSQ[j]);
    37483748          //a = truedist(Lat, 0., Wcentre[j], j); 
    37493749          a = 0.;
    … …  
    37673767        k = cross(i,3);
    37683768        if (GA[g].G[k] == GA[g_bar].G[k]) {
    3769           //b = truedist(Lat, Lat->RealBasisQ[j], Wcentre[j], j);
    3770           b = Lat->RealBasisQ[j];
     3769          //b = truedist(Lat, sqrt(Lat->RealBasisSQ[j]), Wcentre[j], j);
     3770          b = sqrt(Lat->RealBasisSQ[j]);
    37713771          //a = truedist(Lat, 0., Wcentre[j], j); 
    37723772          a = 0.;

  • pcp/src/test.c

    rcc46b0 rf5586e  
    239239        iS = n[2] + N[2]*(n[1] + N[1]*n0);  // mind splitting of x axis due to multiple processes
    240240        i0 = n[2]*NUp[2]+N[2]*NUp[2]*(n[1]*NUp[1]+N[1]*NUp[1]*n0*NUp[0]);
    241         //PsiCR[iS] = ((double)n[cross(index,0)]/(double)N[cross(index,0)]*Lat->RealBasisQ[cross(index,0)] - WCentre[cross(index,0)])*RealPhiR[i0] - ((double)n[cross(index,2)]/(double)N[cross(index,2)]*Lat->RealBasisQ[cross(index,2)] - WCentre[cross(index,2)])*RealPhiR2[i0]; 
     241        //PsiCR[iS] = ((double)n[cross(index,0)]/(double)N[cross(index,0)]*sqrt(Lat->RealBasisSQ[cross(index,0)]) - WCentre[cross(index,0)])*RealPhiR[i0] - ((double)n[cross(index,2)]/(double)N[cross(index,2)]*sqrt(Lat->RealBasisQ[cross(index,2)]) - WCentre[cross(index,2)])*RealPhiR2[i0]; 
    242242        PsiCR[iS] =
    243243            MinImageConv(Lat,x[cross(index,0)],WCentre[cross(index,0)],cross(index,0)) * RealPhiR [i0]

  • pcp/src/wannier.c

    rcc46b0 rf5586e  
    10591059      // calculate Wannier Centre
    10601060      for (j=0;j<NDIM;j++) {
    1061         WannierCentre[i][j] = Lat->RealBasisQ[j]/(2*PI) * GSL_IMAG( gsl_complex_log( gsl_complex_rect(gsl_matrix_get(A[j*2],i,i),gsl_matrix_get(A[j*2+1],i,i))));
     1061        WannierCentre[i][j] = sqrt(Lat->RealBasisSQ[j])/(2*PI) * GSL_IMAG( gsl_complex_log( gsl_complex_rect(gsl_matrix_get(A[j*2],i,i),gsl_matrix_get(A[j*2+1],i,i))));
    10621062        if (WannierCentre[i][j] < 0) // change wrap around of above operator to smooth 0...Lat->RealBasisSQ
    1063           WannierCentre[i][j] = Lat->RealBasisQ[j] + WannierCentre[i][j];
     1063          WannierCentre[i][j] = sqrt(Lat->RealBasisSQ[j]) + WannierCentre[i][j];
    10641064      }
    10651065     
    … …  
    10931093    for (i=0; i < Num; i++) {  // go through all occupied wave functions
    10941094      for (j=0;j<NDIM;j++)
    1095         WannierCentre[i][j] = Lat->RealBasisQ[j]/2.;
     1095        WannierCentre[i][j] = sqrt(Lat->RealBasisSQ[j])/2.;
    10961096    }   
    10971097    break;
    … …  
    11001100    for (i=0;i < Num; i++) {  // go through all wave functions
    11011101      for (j=0;j<NDIM;j++) {
    1102         tmp = WannierCentre[i][j]/Lat->RealBasisQ[j]*(double)N[j];
     1102        tmp = WannierCentre[i][j]/sqrt(Lat->RealBasisSQ[j])*(double)N[j];
    11031103        //fprintf(stderr,"(%i) N[%i]: %i\t tmp %e\t floor %e\t ceil %e\n",P->Par.me, j, N[j], tmp, floor(tmp), ceil(tmp));
    11041104        if (fabs((double)floor(tmp) - tmp) < fabs((double)ceil(tmp) - tmp))
    1105           WannierCentre[i][j] = (double)floor(tmp)/(double)N[j]*Lat->RealBasisQ[j];
     1105          WannierCentre[i][j] = (double)floor(tmp)/(double)N[j]*sqrt(Lat->RealBasisSQ[j]);
    11061106        else
    1107           WannierCentre[i][j] = (double)ceil(tmp)/(double)N[j]*Lat->RealBasisQ[j];
     1107          WannierCentre[i][j] = (double)ceil(tmp)/(double)N[j]*sqrt(Lat->RealBasisSQ[j]);
    11081108      }
    11091109    }
Note: See TracChangeset for help on using the changeset viewer.