Changeset 1d77026 for pcp


Ignore:
Timestamp:
Apr 22, 2008, 11:42:49 AM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
f91abc
Parents:
8786c3
Message:

MinImageConv(): now transform whole vector and does not return single axis value

Due to the switch from hoping for a rectangular simulation box to using Matrix trafo with Real- and ReciBasis, MinImageConv may use these trafos as well and hence may act on the whole vector at once, yielding the minimum image convention in all three directions.
This has impact on numerous functions inside perturbed.c

Location:
pcp/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • pcp/src/perturbed.c

    r8786c3 r1d77026  
    925925  fftw_complex *tempdestRC = (fftw_complex *)Dens0->DensityArray[TempDensity];
    926926  fftw_complex *posfac, *destsnd, *destrcv;
    927   double x[NDIM], fac[NDIM], Wcentre[NDIM];
     927  double x[NDIM], X[NDIM], fac[NDIM], Wcentre[NDIM];
    928928  const int k_normal = Lat->Psi.TypeStartIndex[Occupied] + (wavenr - Lat->Psi.TypeStartIndex[R->CurrentMin]);
    929929  int n[NDIM], n0, g, Index, pos, iS, i0;
    … …  
    998998        //fprintf(stderr,"(%i) R[%i] = (%lg,%lg,%lg)\n",P->Par.me, i0, x[0], x[1], x[2]);
    999999        //else fprintf(stderr,"(%i) WCentre[%i] = %e \n",P->Par.me, index_r, Wcentre[index_r]);
    1000         vector = sawtooth(Lat,MinImageConv(Lat,x[index_r],Wcentre[index_r],index_r),index_r);
     1000        MinImageConv(Lat,x, Wcentre, X);
     1001        vector = sawtooth(Lat,X[index_r],index_r);
    10011002        //vector = 1.;//sin((double)(n[index_r])/(double)((N[index_r]))*2*PI);     
    10021003         PsiCR[iS] = vector * TempPsiR[i0];
    … …  
    12051206  fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity];
    12061207  fftw_real *PsiCR = (fftw_real *) PsiC;
    1207   double x[NDIM], fac[NDIM], Wcentre[NDIM];
     1208  double x[NDIM], X[NDIM], fac[NDIM], Wcentre[NDIM];
    12081209  int n[NDIM], n0, g, Index, iS, i0; //pos,
    12091210  int N[NDIM], NUp[NDIM];
    … …  
    12751276//          }
    12761277       
     1278        MinImageConv(Lat, x, Wcentre, X);
    12771279         PsiCR[iS] = //vector * TempPsiR[i0];
    1278           sawtooth(Lat,MinImageConv(Lat,x[index[0]],Wcentre[index[0]],index[0]),index[0]) * TempPsiR [i0]
    1279          -sawtooth(Lat,MinImageConv(Lat,x[index[2]],Wcentre[index[2]],index[2]),index[2]) * TempPsi2R[i0];
    1280 //          ShiftGaugeOrigin(P,truedist(Lat,x[index[0]],Wcentre[index[0]],index[0]),index[0]) * TempPsiR [i0]
    1281 //         -ShiftGaugeOrigin(P,truedist(Lat,x[index[2]],Wcentre[index[2]],index[2]),index[2]) * TempPsi2R[i0];
     1280          sawtooth(Lat,X[index[0]],index[0]) * TempPsiR [i0]
     1281         -sawtooth(Lat,X[index[2]],index[2]) * TempPsi2R[i0];
     1282//          ShiftGaugeOrigin(P,X[index[0]],index[0]) * TempPsiR [i0]
     1283//         -ShiftGaugeOrigin(P,X[index[2]],index[2]) * TempPsi2R[i0];
    12821284//        PsiCR[iS] = (x[index[0]] - Wcentre[index[0]]) * TempPsiR [i0] - (x[index[2]] - Wcentre[index[2]]) * TempPsi2R[i0];
    12831285      }
    … …  
    13461348  fftw_complex *tempdestRC = (fftw_complex *)Dens0->DensityArray[Temp2Density];
    13471349  fftw_complex *posfac, *destsnd, *destrcv;
    1348   double x[NDIM], fac[NDIM], Wcentre[NDIM];
     1350  double x[NDIM], X[NDIM], fac[NDIM], Wcentre[NDIM];
    13491351  int n[NDIM], n0, g, Index, pos, iS, i0;
    13501352  int N[NDIM], NUp[NDIM];
    … …  
    14041406        iS = n[2] + N[2]*(n[1] + N[1]*n0);  // mind splitting of x axis due to multiple processes
    14051407        i0 = n[2]*NUp[2]+N[2]*NUp[2]*(n[1]*NUp[1]+N[1]*NUp[1]*n0*NUp[0]);
    1406 //         PsiCR[iS] = sawtooth(Lat,truedist(Lat,x[cross(index_pxr,1)],Wcentre[cross(index_pxr,1)],cross(index_pxr,1)),cross(index_pxr,1)) * TempPsiR[i0];
    1407 //         Psi2CR[iS] = sawtooth(Lat,truedist(Lat,x[cross(index_pxr,3)],Wcentre[cross(index_pxr,3)],cross(index_pxr,3)),cross(index_pxr,3)) * TempPsiR[i0];
    1408          PsiCR[iS] = ShiftGaugeOrigin(P,MinImageConv(Lat,x[cross(index_pxr,1)],Wcentre[cross(index_pxr,1)],cross(index_pxr,1)),cross(index_pxr,1)) * TempPsiR[i0];
    1409          Psi2CR[iS] = ShiftGaugeOrigin(P,MinImageConv(Lat,x[cross(index_pxr,3)],Wcentre[cross(index_pxr,3)],cross(index_pxr,3)),cross(index_pxr,3)) * TempPsiR[i0];
     1408//         PsiCR[iS] = sawtooth(Lat,X[cross(index_pxr,1)],cross(index_pxr,1)) * TempPsiR[i0];
     1409//         Psi2CR[iS] = sawtooth(Lat,X[cross(index_pxr,3)],cross(index_pxr,3)) * TempPsiR[i0];
     1410        MinImageConv(Lat,x,Wcentre,X);
     1411        PsiCR[iS] = ShiftGaugeOrigin(P,X[cross(index_pxr,1)],cross(index_pxr,1)) * TempPsiR[i0];
     1412        Psi2CR[iS] = ShiftGaugeOrigin(P,X[cross(index_pxr,3)],cross(index_pxr,3)) * TempPsiR[i0];
    14101413      }
    14111414 
    … …  
    16761679  struct LatticeLevel *LevS = R->LevS;
    16771680  struct Lattice *Lat =&P->Lat;
    1678   double x;
    1679   int n;
     1681  double x[NDIM];
     1682  double n[NDIM];
    16801683  int N[NDIM];
    16811684  N[0] = LevS->Plan0.plan->N[0];
    … …  
    16831686  N[2] = LevS->Plan0.plan->N[2];
    16841687 
    1685   for (n=0;n<N[index];n++) {
    1686     x = (double)n/(double)N[index] * sqrt(Lat->RealBasisSQ[index]);
     1688  n[0] = n[1] = n[2] = 0.;
     1689  for (n[index]=0;n[index]<N[index];n[index]++) {
     1690    n[index] = (double)n[index]/(double)N[index] * sqrt(Lat->RealBasisSQ[index]);
    16871691    //fprintf(stderr,"(%i) x %e\t Axis/2 %e\n",P->Par.me, x, sqrt(Lat->RealBasisSQ[index])/2. );     
    1688     x = MinImageConv(Lat, x, sqrt(Lat->RealBasisSQ[index])/2., index);
    1689     fprintf(stderr,"%e\t%e\n", x, sawtooth(Lat,x,index)); 
    1690   } 
     1692          MinImageConv(Lat, n, Lat->RealBasisCenter, x);
     1693        fprintf(stderr,"%e\t%e\n", n[index], sawtooth(Lat,n[index],index)); 
     1694  }
    16911695}
    16921696
    … …  
    16961700 * \param R[] first vector, NDIM, each must be between 0...L
    16971701 * \param r[] second vector, NDIM, each must be between 0...L
    1698  * \param index of component
    1699  * \return component between  -L/2 ... L/2
     1702 * \param result[] return vector
    17001703 */
    1701 inline double MinImageConv(struct Lattice *Lat, const double R, const double r, const int index)
     1704inline void MinImageConv(struct Lattice *Lat, const double R[NDIM], const double r[NDIM], double *result)
    17021705{
    1703   double axis = sqrt(Lat->RealBasisSQ[index]);
    1704   double result = 0.;
    1705 
    1706   if (fabs(result = R - r + axis) < axis/2.) { }
    1707   else if (fabs(result = R - r) <= axis/2.) { }
    1708    else if (fabs(result = R - r - axis) < axis/2.) { }
    1709     else Error(SomeError, "MinImageConv: None of the three cases applied!");
    1710   return (result);
     1706  //double axis = Lat->RealBasisQ[index];
     1707  double x[NDIM], X[NDIM], Result[NDIM];
     1708  int i;
     1709 
     1710  for(i=0;i<NDIM;i++)
     1711    result[i] = x[i] = x[i] = 0.;
     1712  //fprintf(stderr, "R = (%lg, %lg, %lg), r = (%lg, %lg, %lg)\n", R[0], R[1], R[2], r[0], r[1], r[2]);
     1713  RMat33Vec3(X, Lat->ReciBasis, R); // transform both to [0,1]^3
     1714  RMat33Vec3(x, Lat->ReciBasis, r);
     1715  //fprintf(stderr, "X = (%lg, %lg, %lg), x = (%lg, %lg, %lg)\n", X[0], X[1], X[2], x[0], x[1], x[2]);
     1716  for(i=0;i<NDIM;i++) {
     1717//    if (fabs(X[i]) > 1.)
     1718//      fprintf(stderr,"X[%i] > 1. : %lg!\n", i, X[i]);
     1719//    if (fabs(x[i]) > 1.)
     1720//      fprintf(stderr,"x[%i] > 1. : %lg!\n", i, x[i]);
     1721    if (fabs(Result[i] = X[i] - x[i] + 2.*PI) < PI) { }
     1722    else if (fabs(Result[i] = X[i] - x[i]) <= PI) { }
     1723     else if (fabs(Result[i] = X[i] - x[i] - 2.*PI) < PI) { }
     1724      else Error(SomeError, "MinImageConv: None of the three cases applied!");
     1725  }
     1726  for(i=0;i<NDIM;i++) // ReciBasis is not true inverse, but times 2.*PI
     1727    Result[i] /= 2.*PI;
     1728  RMat33Vec3(result, Lat->RealBasis, Result);
    17111729}
    17121730
    … …  
    20852103  fftw_real *Psip1R;
    20862104  fftw_real *tempArray; // intendedly the same
    2087   double r_bar[NDIM], x[NDIM], fac[NDIM];
     2105  double r_bar[NDIM], x[NDIM], X[NDIM], fac[NDIM];
    20882106  double Current;//, current;
    20892107  const double UnitsFactor = 1.; ///LevS->MaxN;  // 1/N (from ff-backtransform)
    … …  
    22492267                fac[2] = (double)n[2]/(double)N[2];
    22502268                RMat33Vec3(x, Lat->RealBasis, fac); // relative coordinate times basis matrix gives absolute ones
     2269                MinImageConv(Lat, x, Psi->AddData[k].WannierCentre, X);
    22512270                for (i=0;i<NDIM;i++) // build gauge-translated r_bar evaluation point
    2252                   r_bar[i] =
    2253                     sawtooth(Lat,MinImageConv(Lat, x[i], Psi->AddData[k].WannierCentre[i], i),i);
    2254 //                    ShiftGaugeOrigin(P,truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i),i);
    2255                     //truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i);
     2271                  r_bar[i] = sawtooth(Lat,X[i],i);
     2272//                    ShiftGaugeOrigin(P,X[i],i);
     2273                    //X[i];
    22562274                Current = Psip0R[i0] * (r_bar[cross_lookup_1[0]] * Psi1R[i0]);
    22572275                Current += (Psi0R[i0] * r_bar[cross_lookup_1[0]] * Psip1R[i0]);
    … …  
    25012519  int mem_avail, MEM_avail;
    25022520  double Current;
     2521  double X[NDIM];
    25032522  const double UnitsFactor = 1.;
    25042523  int cross_lookup[4];
    … …  
    26312650                  for (dex=0;dex<4;dex++)
    26322651                    cross_lookup[dex] = cross(in,dex);
     2652                  MinImageConv(Lat,Psi->AddData[LOnePsiA->MyLocalNo].WannierCentre, Psi->AddData[LOnePsiB->MyLocalNo].WannierCentre,X);
    26332653                  for(s=0;s<LevS->MaxG;s++) {
    26342654                    //if (x_l != x_l_bak || s<0 || s>LevS->MaxG) Error(SomeError,"FillDeltaCurrentDensity: x_l[] corrupted");
    2635                     factor =
    2636                       (MinImageConv(Lat,Psi->AddData[LOnePsiA->MyLocalNo].WannierCentre[cross_lookup[0]],
    2637                                         Psi->AddData[LOnePsiB->MyLocalNo].WannierCentre[cross_lookup[0]],cross_lookup[0]) * LevS->GArray[s].G[cross_lookup[1]] -
    2638                        MinImageConv(Lat,Psi->AddData[LOnePsiA->MyLocalNo].WannierCentre[cross_lookup[2]],
    2639                                         Psi->AddData[LOnePsiB->MyLocalNo].WannierCentre[cross_lookup[2]],cross_lookup[2]) * LevS->GArray[s].G[cross_lookup[3]]);
     2655                    factor = (X[cross_lookup[0]] * LevS->GArray[s].G[cross_lookup[1]] - X[cross_lookup[2]] * LevS->GArray[s].G[cross_lookup[3]]);
    26402656                    x_l[s].re = factor * (-LPsiDatB[s].im); // switched due to factorization with "-i G"
    26412657                    x_l[s].im = factor * (LPsiDatB[s].re);
    … …  
    28772893  N[1] = Lev0->Plan0.plan->N[1];
    28782894  N[2] = Lev0->Plan0.plan->N[2];
    2879   double chi[NDIM*NDIM],Chi[NDIM*NDIM], x[NDIM], fac[NDIM];
     2895  double chi[NDIM*NDIM],Chi[NDIM*NDIM], x[NDIM], X[NDIM], fac[NDIM];
    28802896  const double discrete_factor = Lat->Volume/Lev0->MaxN;
    28812897  const int myPE =  P->Par.me_comm_ST_Psi;
    … …  
    29212937            RMat33Vec3(x, Lat->RealBasis, fac);
    29222938            i0 = n[2]+N[2]*(n[1]+N[1]*(n0));  // the index of current density must match LocalSizeR!
    2923             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.
    2924             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.
     2939            MinImageConv(Lat,x, Lat->RealBasisCenter, X);
     2940            chi[in+dex*NDIM] += X[cross_lookup[0]] * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]; // x[cross(in,0)], Lat->RealBasisCenter[cross_lookup[0]]
     2941            chi[in+dex*NDIM] -= X[cross_lookup[2]] * CurrentDensity[dex*NDIM+cross_lookup[3]][i0]; // x[cross(in,2)], Lat->RealBasisCenter[cross_lookup[2]]
    29252942//            if (in == dex) field[in][i0] =
    29262943//                truedist(Lat,x[cross_lookup[0]], sqrt(Lat->RealBasisSQ[c[0]])/2.,cross_lookup[0]) * CurrentDensity[dex*NDIM+cross_lookup[1]][i0]
    … …  
    30223039  fftw_real *sigma_real[NDIM_NDIM];
    30233040  double sigma,Sigma;
    3024   int it, ion, in, dex, g, n[2][NDIM], Index, i;
     3041  double x[2][NDIM];
     3042  int it, ion, in, dex, g, Index, i;
    30253043  //const double FFTfactor = 1.;///Lev0->MaxN;
     3044  int n[2][NDIM];
    30263045  double eta, delta_sigma, S, A, iso, tmp;
    30273046  FILE *SigmaFile;
    30283047  char suffixsigma[255];
    3029   double x[2][NDIM];
    30303048  const int myPE = P->Par.me_comm_ST_Psi;
    30313049  int N[NDIM];
    … …  
    32723290  fftw_real *CurrentDensity[NDIM*NDIM];
    32733291  int it, ion, in, dex, i0, n[NDIM], n0, i;//, *NUp;
    3274   double x,y,z;
    3275   double dist;
    3276   double r[NDIM], fac[NDIM];
     3292  double r[NDIM], fac[NDIM], X[NDIM];
    32773293  const double discrete_factor = Lat->Volume/Lev0->MaxN;
    32783294  double eta, delta_sigma, S, A, iso;
    … …  
    33153331                fac[2] = (double)n[2]/(double)N[2];
    33163332                RMat33Vec3(r, Lat->RealBasis, fac);
     3333                MinImageConv(Lat,r, &(I->I[it].R[NDIM*ion]),X);
    33173334                i0 = n[2]+N[2]*(n[1]+N[1]*(n0));  // the index of current density must match LocalSizeR!
    3318                 x = MinImageConv(Lat,r[cross(in,0)], I->I[it].R[NDIM*ion+cross(in,0)],cross(in,0));
    3319                 y = MinImageConv(Lat,r[cross(in,2)], I->I[it].R[NDIM*ion+cross(in,2)],cross(in,2));
    3320                 z = MinImageConv(Lat,r[in], I->I[it].R[NDIM*ion+in],in);  // "in" always is missing third component in cross product
    3321                 dist = pow(x*x + y*y + z*z, 3./2.);
    3322                 if ((dist < pow(2.,3.)) && (dist > MYEPSILON)) sigma[in+dex*NDIM] += (x * CurrentDensity[dex*NDIM+cross(in,1)][i0] - y * CurrentDensity[dex*NDIM+cross(in,3)][i0])/dist;
    3323                 //if (it == 0 && ion == 0) fprintf(stderr,"(%i) sigma[%i][%i] += (%e * %e - %e * %e)/%e = %e\n", P->Par.me, in, dex, x,CurrentDensity[dex*NDIM+cross(in,1)][i0],y,CurrentDensity[dex*NDIM+cross(in,3)][i0],dist,sigma[in+dex*NDIM]);
     3335                //z = MinImageConv(Lat,r, I->I[it].R[NDIM*ion],in);  // "in" always is missing third component in cross product
     3336                sigma[in+dex*NDIM] += (X[cross(in,0)] * CurrentDensity[dex*NDIM+cross(in,1)][i0] - X[cross(in,2)] * CurrentDensity[dex*NDIM+cross(in,3)][i0]);
     3337                //if (it == 0 && ion == 0) fprintf(stderr,"(%i) moment[%i][%i] += (%e * %e - %e * %e) = %e\n", P->Par.me, in, dex, x,CurrentDensity[dex*NDIM+cross(in,1)][i0],y,CurrentDensity[dex*NDIM+cross(in,3)][i0],moment[in+dex*NDIM]);
    33243338              }
    33253339          sigma[in+dex*NDIM] *= -mu0*discrete_factor/(4.*PI); // due to summation instead of integration

  • pcp/src/perturbed.h

    r8786c3 r1d77026  
    4444inline double ShiftGaugeOrigin(struct Problem *P, double r, const int index);
    4545inline double sawtooth(struct Lattice *Lat, const double L, const int index);
    46 inline double MinImageConv(struct Lattice *Lat, const double R, const double r, const int index);
     46inline void MinImageConv(struct Lattice *Lat, const double R[NDIM], const double r[NDIM], double *result);
    4747void test_fft_symmetry(struct Problem *P, const int l);
    4848void test_rxp(struct Problem *P);
Note: See TracChangeset for help on using the changeset viewer.