/*
    DFT++ is a density functional package developed by the research group
    of Professor Tomas Arias

    Copyright 1996-2003 Sohrab Ismail-Beigi

    This file is part of DFT++.

    DFT++ is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    DFT++ is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DFT++; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Please see the file CREDITS for a list of authors.

    For academic users, we request that publications using results obtained with
    this software reference

    "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
    and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).

    and, if using the wavelet basis, further reference

    "Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
    T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).

    and 

    "Robust ab initio calculation of condensed matter: transparent convergence through
    semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
    Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).

    For your convenience, preprints of the above articles may be obtained from
    http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/

#include "header.h"

/*
 * check_symm_fftbox
 *
 * Make sure that the FFT-box doesn't destroy the intrinsic symmetries of 
 * the system.
 *
 * And calculate the symmetry matrix for charge density mesh.
 */
int check_basis_symm_compatibility(Symmetries &symm, const PW_Basis &basis)
{
  int N[3] = { basis.basis_spec->Nx,
               basis.basis_spec->Ny,
               basis.basis_spec->Nz };
  matrix3 tmp, Nmat(basis.basis_spec->Nx,
                    basis.basis_spec->Ny,
                    basis.basis_spec->Nz),
    Nmat_inv(1.0/basis.basis_spec->Nx,
             1.0/basis.basis_spec->Ny,
             1.0/basis.basis_spec->Nz);
  int nrot = symm.nrot, irot, i1, i2, error = FALSE;
  
  for (irot = 0; irot < nrot; irot++) {
    tmp = symm.sym[irot] * Nmat;
    for (i1 = 0; i1 < 3; i1++) {
      for (i2 = 0; i2 < 3; i2++) {
        if (fmod(tmp.m[i1][i2], N[i1]) > MIN_SYMM_TOL) {
          if (!error) 
            dft_log(">>>> FFT box not comensurate with symmetries:\n");
          symm.sym[irot].print(dft_global_log,"%4.0f ");
          dft_log(">>>> \n");
          error = TRUE;
        }
      }
    }
    symm.n_sym[irot] = Nmat * symm.sym[irot] * Nmat_inv;
  }
  if (error) {
    dft_log(DFT_SILENCE,">>>> \t[ %d %d %d ]\n",N[0], N[1], N[2]);
    dft_log(DFT_SILENCE,">>>> You have two options:\n\t1. disable symmetry\n");
    dft_log(DFT_SILENCE,"\t2. manually input FFT box size\n");
    die("");
  }

  return (! error);

}

/*
 * Symmetrize the charge density for planewaves
 */
int symmetrize_n(RealSpaceScalarFieldColumn &n, Symmetries *symm)
{
  Basis * basis = n.basis;
  // Z is the outer index, Y the next, X the inner index.
  vector3 r, rnew;
  int irot, i;
  int index, index2, N[3], ind[48];
  real one_over_nrot, n_sum;

  if (symm->nrot <= 1)  // no need to do any symmetrization
    return 0;

  one_over_nrot = 1.0/symm->nrot;
  for (i = 0; i < n.length; i++)
    symm->done[i] = FALSE;

  N[0] = basis->basis_spec->Nx;
  N[1] = basis->basis_spec->Ny;
  N[2] = basis->basis_spec->Nz;
  
  index = 0;
  
  for (r.v[0]=0; r.v[0] < N[0]; r.v[0]+=1) {
    for (r.v[1]=0; r.v[1] < N[1]; r.v[1]+=1) {
      for (r.v[2]=0; r.v[2] < N[2]; r.v[2]+=1, index++) {
        // for all points in the box
        
        // if not yet related by symmetry to a previous point
        if ( !symm->done[index] ) {
          
          n_sum = REAL( n.data.d[index] );   // first symmetry is identity.
          ind[0] = index;
          // for all other symmetries
          for (irot = 1; irot < symm->nrot; irot++) {
            // get the point related by symmetry.
            rnew = symm->n_sym[irot] * r;
            
            for (i=0;i<3;i++) {	    // project back into range.
              rnew.v[i] = rint(rnew.v[i]);
              rnew.v[i] = fmod(rnew.v[i], N[i]);
              if (rnew.v[i] < 0.0)
                rnew.v[i] += N[i];
            }
            index2 = (int)(rnew.v[2] + N[2]*(rnew.v[1] + N[1]*rnew.v[0]));
            
            n_sum += REAL(n.data.d[ index2 ]); // add the charge density.
            
            ind[irot] = index2; // remember the index.
          }
          
          n_sum *= one_over_nrot;  // get averaged charge density.
          
          // put symmetrized value in the star.
          for (irot = 0; irot < symm->nrot; irot++) {
            n.data.d[ ind[irot] ] = n_sum;
            symm->done[ ind[irot] ] = TRUE;
          }
        }
      }
    }
  }

  return 1;
}

int symmetrize_n_dag(RealSpaceScalarFieldColumn &n, Symmetries *symm)
{
  return 0;
}