fields.c

麻省理工的计算光子晶体的程序

#endif
     }

     mpi_allreduce_1(&energy_sum, real, SCALAR_MPI_TYPE,
		     MPI_SUM, MPI_COMM_WORLD);
     mpi_allreduce(comp_sum2, comp_sum, 6, real, SCALAR_MPI_TYPE,
                   MPI_SUM, MPI_COMM_WORLD);

     mpi_one_printf("%c-energy-components:, %d, %d",
	    curfield_type, kpoint_index, curfield_band);
     for (i = 0; i < 6; ++i) {
	  comp_sum[i] /= (energy_sum == 0 ? 1 : energy_sum);
	  if (i % 2 == 1)
	       mpi_one_printf(", %g", comp_sum[i] + comp_sum[i-1]);
     }
     mpi_one_printf("\n");


     /* remember that we now have energy density; denoted by capital D/H */
     curfield_type = toupper(curfield_type);

     /* The return value is a list of 7 items: the total energy,
	followed by the 6 elements of the comp_sum array (the fraction
	of the energy in the real/imag. parts of each field component). */

     retval.num_items = 7;
     CHK_MALLOC(retval.items, number, retval.num_items);

     retval.items[0] = energy_sum * Vol / H.N;

     for (i = 0; i < 6; ++i)
	  retval.items[i+1] = comp_sum[i];

     return retval;
}

/**************************************************************************/

/* Fix the phase of the current field (e/h/d) to a canonical value.
   Also changes the phase of the corresponding eigenvector by the
   same amount, so that future calculations will have a consistent
   phase.

   The following procedure is used, derived from a suggestion by Doug
   Allan of Corning: First, choose the phase to maximize the sum of
   the squares of the real parts of the components.  This doesn't fix
   the overall sign, though.  That is done (after incorporating the
   above phase) by: (1) find the largest absolute value of the real
   part, (2) find the point with the greatest spatial array index that
   has |real part| at least half of the largest value, and (3) make
   that point positive.

   In the case of inversion symmetry, on the other hand, the overall phase
   is already fixed, to within a sign, by the choice to make the Fourier
   transform purely real.  So, in that case we simply pick a sign, in
   a manner similar to (2) and (3) above. */
void fix_field_phase(void)
{
     int i, N;
     real sq_sum2[2] = {0,0}, sq_sum[2], maxabs = 0.0;
     int maxabs_index = 0, maxabs_sign = 1;
     double theta;
     scalar phase;

     if (!curfield || !strchr("dhecv", curfield_type)) {
          mpi_one_fprintf(stderr, "The D/H/E field must be loaded first.\n");
          return;
     }
     N = mdata->fft_output_size * 3;

#ifdef SCALAR_COMPLEX
     /* Compute the phase that maximizes the sum of the squares of
	the real parts of the components.  Equivalently, maximize
	the real part of the sum of the squares. */
     for (i = 0; i < N; ++i) {
	  real a,b;
	  a = curfield[i].re; b = curfield[i].im;
	  sq_sum2[0] += a*a - b*b;
	  sq_sum2[1] += 2*a*b;
     }
     mpi_allreduce(sq_sum2, sq_sum, 2, real, SCALAR_MPI_TYPE,
                   MPI_SUM, MPI_COMM_WORLD);
     /* compute the phase = exp(i*theta) maximizing the real part of
	the sum of the squares.  i.e., maximize:
	    cos(2*theta)*sq_sum[0] - sin(2*theta)*sq_sum[1] */
     theta = 0.5 * atan2(-sq_sum[1], sq_sum[0]);
     phase.re = cos(theta);
     phase.im = sin(theta);
#else /* ! SCALAR_COMPLEX */
     phase = 1;
#endif /* ! SCALAR_COMPLEX */

     /* Next, fix the overall sign.  We do this by first computing the
	maximum |real part| of the jmax component (after multiplying
	by phase), and then finding the last spatial index at which
	|real part| is at least half of this value.  The sign is then
	chosen to make the real part positive at that point. 

        (Note that we can't just make the point of maximum |real part|
         positive, as that would be ambiguous in the common case of an
         oscillating field within the unit cell.)

        In the case of inversion symmetry (!SCALAR_COMPLEX), we work with
        (real part - imag part) instead of (real part), to insure that we
        have something that is nonzero somewhere. */

     for (i = 0; i < N; ++i) {
#ifdef SCALAR_COMPLEX
	  real r = fabs(curfield[i].re * phase.re - curfield[i].im * phase.im);
#else
	  real r = fabs(curfield[i].re - curfield[i].im);
#endif
	  if (r > maxabs)
	       maxabs = r;
     }
     mpi_allreduce_1(&maxabs, real, SCALAR_MPI_TYPE,
		     MPI_MAX, MPI_COMM_WORLD);
     for (i = N - 1; i >= 0; --i) {
#ifdef SCALAR_COMPLEX
	  real r = curfield[i].re * phase.re - curfield[i].im * phase.im;
#else
	  real r = curfield[i].re - curfield[i].im;
#endif
	  if (fabs(r) >= 0.5 * maxabs) {
	       maxabs_index = i;
	       maxabs_sign = r < 0 ? -1 : 1;
	       break;
	  }
     }
     if (i >= 0)  /* convert index to global index in distributed array: */
	  maxabs_index += mdata->local_y_start * mdata->nx * mdata->nz;
     {
	  /* compute maximum index and corresponding sign over all the 
	     processors, using the MPI_MAXLOC reduction operation: */
	  struct twoint_struct {int i; int s;} x;
	  x.i = maxabs_index; x.s = maxabs_sign;
	  mpi_allreduce_1(&x, struct twoint_struct, MPI_2INT,
			  MPI_MAXLOC, MPI_COMM_WORLD);
	  maxabs_index = x.i; maxabs_sign = x.s;
     }
     ASSIGN_SCALAR(phase,
		   SCALAR_RE(phase)*maxabs_sign, SCALAR_IM(phase)*maxabs_sign);

     mpi_one_printf("Fixing %c-field (band %d) phase by %g + %gi; "
		    "max ampl. = %g\n", curfield_type, curfield_band,
		    SCALAR_RE(phase), SCALAR_IM(phase), maxabs);

     /* Now, multiply everything by this phase, *including* the
	stored "raw" eigenvector in H, so that any future fields
	that we compute will have a consistent phase: */
     for (i = 0; i < N; ++i) {
	  real a,b;
	  a = curfield[i].re; b = curfield[i].im;
	  curfield[i].re = a*SCALAR_RE(phase) - b*SCALAR_IM(phase);
	  curfield[i].im = a*SCALAR_IM(phase) + b*SCALAR_RE(phase);
     }
     for (i = 0; i < H.n; ++i) {
          ASSIGN_MULT(H.data[i*H.p + curfield_band - 1], 
		      H.data[i*H.p + curfield_band - 1], phase);
     }
}

/**************************************************************************/

/* Functions to return epsilon, fields, energies, etcetera, at a specified
   point, linearly interpolating if necessary. */

static real get_val(int ix, int iy, int iz,
		    int nx, int ny, int nz, int last_dim_size,
		    real *data, int stride, int conjugate)
{
#ifndef SCALAR_COMPLEX
     {
	  int nlast = last_dim_size / 2;
	  if ((nz > 1 ? iz : (ny > 1 ? iy : ix)) >= nlast) {
	       ix = ix ? nx - ix : ix;
	       iy = iy ? ny - iy : iy;
	       iz = iz ? nz - iz : iz;
	       conjugate = conjugate ? 1 : 0;
	  }
	  else
	       conjugate = 0;
	  if (nz > 1) nz = nlast; else if (ny > 1) ny = nlast; else nx = nlast;
     }
#else
     conjugate = 0;
#endif

#ifdef HAVE_MPI
     CHECK(0, "get-*-point not yet implemented for MPI!");
#else
     if (conjugate)
	  return -data[(((ix * ny) + iy) * nz + iz) * stride];
     else
	  return data[(((ix * ny) + iy) * nz + iz) * stride];
#endif
}

static real interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size,
		       real *data, int stride, int conjugate)
{
     double ipart;
     real rx, ry, rz, dx, dy, dz;
     int x, y, z, x2, y2, z2;

     rx = modf(p.x/geometry_lattice.size.x + 0.5, &ipart); if (rx < 0) rx += 1;
     ry = modf(p.y/geometry_lattice.size.y + 0.5, &ipart); if (ry < 0) ry += 1;
     rz = modf(p.z/geometry_lattice.size.z + 0.5, &ipart); if (rz < 0) rz += 1;

     /* get the point corresponding to r in the grid: */
     x = rx * nx;
     y = ry * ny;
     z = rz * nz;

     /* get the difference between (x,y,z) and the actual point */
     dx = rx * nx - x;
     dy = ry * ny - y;
     dz = rz * nz - z;

     /* get the other closest point in the grid, with periodic boundaries: */
     x2 = (nx + (dx >= 0.0 ? x + 1 : x - 1)) % nx;
     y2 = (ny + (dy >= 0.0 ? y + 1 : y - 1)) % ny;
     z2 = (nz + (dz >= 0.0 ? z + 1 : z - 1)) % nz;

     /* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
     dx = fabs(dx);
     dy = fabs(dy);
     dz = fabs(dz);

#define D(x,y,z) (get_val(x,y,z,nx,ny,nz,last_dim_size, data,stride,conjugate))

     return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
             (D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
            ((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
             (D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);

#undef D
}

static scalar_complex interp_cval(vector3 p, 
				  int nx, int ny, int nz, int last_dim_size,
				  real *data, int stride)
{
     scalar_complex cval;
     cval.re = interp_val(p, nx,ny,nz,last_dim_size, data, stride, 0);
     cval.im = interp_val(p, nx,ny,nz,last_dim_size,data + 1, stride, 1);
     return cval;
}

#define f_interp_val(p,f,data,stride,conj) interp_val(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride,conj)
#define f_interp_cval(p,f,data,stride) interp_cval(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride)

static symmetric_matrix interp_eps_inv(vector3 p)
{
     int stride = sizeof(symmetric_matrix) / sizeof(real);
     symmetric_matrix eps_inv;

     eps_inv.m00 = f_interp_val(p, mdata, &mdata->eps_inv->m00, stride, 0);
     eps_inv.m11 = f_interp_val(p, mdata, &mdata->eps_inv->m11, stride, 0);
     eps_inv.m22 = f_interp_val(p, mdata, &mdata->eps_inv->m22, stride, 0);
#ifdef WITH_HERMITIAN_EPSILON
     eps_inv.m01 = f_interp_cval(p, mdata, &mdata->eps_inv->m01.re, stride);
     eps_inv.m02 = f_interp_cval(p, mdata, &mdata->eps_inv->m02.re, stride);
     eps_inv.m12 = f_interp_cval(p, mdata, &mdata->eps_inv->m12.re, stride);
#else
     eps_inv.m01 = f_interp_val(p, mdata, &mdata->eps_inv->m01, stride, 0);
     eps_inv.m02 = f_interp_val(p, mdata, &mdata->eps_inv->m02, stride, 0);
     eps_inv.m12 = f_interp_val(p, mdata, &mdata->eps_inv->m12, stride, 0);
#endif
     return eps_inv;
}

number get_epsilon_point(vector3 p)
{
     symmetric_matrix eps_inv;
     eps_inv = interp_eps_inv(p);
     return mean_epsilon(&eps_inv);
}

cmatrix3x3 get_epsilon_inverse_tensor_point(vector3 p)
{
     symmetric_matrix eps_inv;
     eps_inv = interp_eps_inv(p);

#ifdef WITH_HERMITIAN_EPSILON
     return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
				      cscalar2cnumber(eps_inv.m01),
				      cscalar2cnumber(eps_inv.m02),
				      cscalar2cnumber(eps_inv.m12));
#else
     return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
				      make_cnumber(eps_inv.m01,0),
				      make_cnumber(eps_inv.m02,0),
				      make_cnumber(eps_inv.m12,0));
#endif
}

number get_energy_point(vector3 p)
{
     CHECK(curfield && strchr("DHR", curfield_type),
	   "compute-field-energy must be called before get-energy-point");
     return f_interp_val(p, mdata, (real *) curfield, 1, 0);
}

cvector3 get_bloch_field_point(vector3 p)
{
     scalar_complex field[3];
     cvector3 F;

     CHECK(curfield && strchr("dhecv", curfield_type),
	   "field must be must be loaded before get-*field*-point");
     field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
     field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
     field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);
     F.x = cscalar2cnumber(field[0]);
     F.y = cscalar2cnumber(field[1]);
     F.z = cscalar2cnumber(field[2]);
     return F;
}

cvector3 get_field_point(vector3 p)
{
     scalar_complex field[3], phase;
     cvector3 F;

     CHECK(curfield && strchr("dhecv", curfield_type),
	   "field must be must be loaded before get-*field*-point");
     field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
     field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
     field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);

     if (curfield_type != 'v') {
	  double phase_phi = TWOPI * 
	       (cur_kvector.x * (p.x/geometry_lattice.size.x+0.5) +
		cur_kvector.y * (p.y/geometry_lattice.size.y+0.5) +
		cur_kvector.z * (p.z/geometry_lattice.size.z+0.5));
	  CASSIGN_SCALAR(phase, cos(phase_phi), sin(phase_phi));
	  CASSIGN_MULT(field[0], field[0], phase);
	  CASSIGN_MULT(field[1], field[1], phase);
	  CASSIGN_MULT(field[2], field[2], phase);
     }

     F.x = cscalar2cnumber(field[0]);
     F.y = cscalar2cnumber(field[1]);
     F.z = cscalar2cnumber(field[2]);
     return F;
}

number rscalar_field_get_point(SCM fo, vector3 p)
{
     field_smob *f = assert_field_smob(fo);
     CHECK(f->type == RSCALAR_FIELD_SMOB, 
	   "invalid argument to rscalar-field-get-point");
     return f_interp_val(p, f, f->f.rs, 1, 0);
}

cvector3 cvector_field_get_point_bloch(SCM fo, vector3 p)
{
     scalar_complex field[3];
     cvector3 F;
     field_smob *f = assert_field_smob(fo);
     CHECK(f->type == CVECTOR_FIELD_SMOB, 
	   "invalid argument to cvector-field-get-point");
     field[0] = f_interp_cval(p, f, &f->f.cv[0].re, 6);
     field[1] = f_interp_cval(p, f, &f->f.cv[1].re, 6);
     field[2] = f_interp_cval(p, f, &f->f.cv[2].re, 6);
     F.x = cscalar2cnumber(field[0]);
     F.y = cscalar2cnumber(field[1]);
     F.z = cscalar2cnumber(field[2]);
     return F;
}

cvector3 cvector_field_get_point(SCM fo, vector3 p)
{
     scalar_complex field[3];
     cvector3 F;
     field_smob *f = assert_field_smob(fo);
     CHECK(f->type == CVECTOR_FIELD_SMOB, 
	   "invalid argument to cvector-field-get-point");

     field[0] = f_interp_cval(p, f, &f->f.cv[0].re, 6);
     field[1] = f_interp_cval(p, f, &f->f.cv[1].re, 6);
     field[2] = f_interp_cval(p, f, &f->f.cv[2].re, 6);
     
     if (f->type_char != 'v') { /* v fields have no kvector */
	  scalar_complex phase;
	  double phase_phi = TWOPI * 
	       (cur_kvector.x * (p.x/geometry_lattice.size.x+0.5) +
		cur_kvector.y * (p.y/geometry_lattice.size.y+0.5) +
		cur_kvector.z * (p.z/geometry_lattice.size.z+0.5));
	  CASSIGN_SCALAR(phase, cos(phase_phi), sin(phase_phi));
	  CASSIGN_MULT(field[0], field[0], phase);
	  CASSIGN_MULT(field[1], field[1], phase);