Changeset 0da6d5 for pcp

Timestamp:

Apr 22, 2008, 11:11:04 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

d2f1b1

Parents:

b0225a

Message:

FillCurrentDensity(): N is now ptr not array, commented-out Factor part included, verbosity changes

N that holds number of nodes per axis now just points to LatticeLevel structure instead of copying values.
Factor part was meant to balance P- against RxP-Perturbation, but proved futile so far.
some if (P->Call.out[ReadOut]) added

File:

• pcp/src/perturbed.c (modified) (6 diffs)

pcp/src/perturbed.c

@@ -2091 +2091 @@
   int k, j, i0;
   int n[NDIM], n0;
-  int N[NDIM];
-  N[0] = Lev0->Plan0.plan->N[0];
-  N[1] = Lev0->Plan0.plan->N[1];
-  N[2] = Lev0->Plan0.plan->N[2];
+  int *N;
+  N = Lev0->Plan0.plan->N;
   const int N0 = Lev0->Plan0.plan->local_nx;
   //int ActNum;
@@ -2101 +2099 @@
   MPI_Status status;
   int cross_lookup_1[4], cross_lookup_3[4], l_1 = 0, l_3 = 0;
+  double Factor;//, factor;
 
   //fprintf(stderr,"(%i) FactoR %e\n", P->Par.me, R->FactorDensityR);
@@ -2175 +2174 @@
   int wished = -1;
   FILE *file =  fopen(P->Call.MainParameterFile,"r");
-  if (!ParseForParameter(1,file,"Orbital",0,1,1,int_type,&wished, 1, optional)) {
-    fprintf(stderr,"Desired Orbital missing!\n");
+  if (!ParseForParameter(0,file,"Orbital",0,1,1,int_type,&wished, 1, optional)) {
+     if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital missing, using: All!\n");
     wished = -1;
   } else if (wished != -1) {
-    fprintf(stderr,"Desired Orbital is: %i.\n", wished);
+     if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital is: %i.\n", wished);
   } else { 
-    fprintf(stderr,"Desired Orbital is: All.\n");
+     if (P->Call.out[ReadOut]) fprintf(stderr,"Desired Orbital is: All.\n");
   }
   fclose(file);
@@ -2218 +2217 @@
         case Perturbed_P1:
         case Perturbed_P2:
+/*          // evaluate factor to compensate r x normalized phi(r) against normalized phi(rxp)
+          factor = 0.;
+          for (n0=0;n0<N0;n0++)  // only local points on x axis
+            for (n[1]=0;n[1]<N[1];n[1]++)
+              for (n[2]=0;n[2]<N[2];n[2]++) {
+                i0 = n[2]+N[2]*(n[1]+N[1]*n0);
+                n[0]=n0 + N0*myPE; // global relative coordinate: due to partitoning of x-axis in PEPGamma>1 case
+                fac[0] = (double)n[0]/(double)N[0];
+                fac[1] = (double)n[1]/(double)N[1];
+                fac[2] = (double)n[2]/(double)N[2];
+                RMat33Vec3(x, Lat->RealBasis, fac); // relative coordinate times basis matrix gives absolute ones
+                MinImageConv(Lat, x, Psi->AddData[k].WannierCentre, X)
+                for (i=0;i<NDIM;i++) // build gauge-translated r_bar evaluation point
+                  r_bar[i] = sawtooth(Lat,X,i);
+//                    ShiftGaugeOrigin(P,truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i),i);
+                    //truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i);
+                factor += Psi1R[i0] * (r_bar[cross_lookup_1[0]] * Psi1R[i0]);
+              }
+          MPI_Allreduce (&factor, &Factor, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
+          Factor *= R->FactorDensityR;  // discrete integration constant
+          fprintf(stderr,"(%i) normalization factor of Phi^(RxP%i)_{%i} is %lg\n", P->Par.me, type, k, Factor);
+          Factor =  1./sqrt(fabs(Factor)); //Factor/fabs(Factor) */
+          Factor = 1.;
           for (n0=0;n0<N0;n0++)  // only local points on x axis
             for (n[1]=0;n[1]<N[1];n[1]++)
@@ -2232 +2254 @@
 //                    ShiftGaugeOrigin(P,truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i),i);
                     //truedist(Lat, x[i], Psi->AddData[k].WannierCentre[i], i);
-                
                 Current = Psip0R[i0] * (r_bar[cross_lookup_1[0]] * Psi1R[i0]);
                 Current += (Psi0R[i0] * r_bar[cross_lookup_1[0]] * Psip1R[i0]);
@@ -2238 +2259 @@
                 ////if (CurrentDensity[index+j*NDIM] != (fftw_real *) Dens0->DensityArray[CurrentDensity0 + index+j*NDIM] || i0<0 || i0>=Dens0->LocalSizeR || (index+j*NDIM)<0 || (index+j*NDIM)>=NDIM*NDIM) Error(SomeError,"FillCurrentDensity: CurrentDensity[] corrupted");
                 CurrentDensity[index+l_1*NDIM][i0] -= Current; // note: sign of cross product resides in Current itself (here: plus)
-                
                 Current = - Psip0R[i0] * (r_bar[cross_lookup_3[2]] * Psi1R[i0]);
                 Current += - (Psi0R[i0] * r_bar[cross_lookup_3[2]] * Psip1R[i0]);