/*
    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"

// constructor
Ioninfo::Ioninfo()
{
  nspecies = 0;
  species = NULL;
  read_vloc_flag = FALSE;
}

Ioninfo::~Ioninfo()
{
  if (nspecies>0)
    {
      for (int sp=0; sp < nspecies; sp++)
	species[sp].Speciesinfo::~Speciesinfo();
      myfree(species);
    }
}

//
// Do the setup for the ionic information (read psps, alloc space, blah blah)
// 
void
Ioninfo::setup(Everything &e)
{
  dft_log("----- Ioninfo::setup() -----\n");

  // Call the species setup routines
  for (int sp=0; sp < e.ioninfo.nspecies; sp++)
    e.ioninfo.species[sp].setup(e);

  // Check atomic positions for problems
  check_pos();

  // Print out the internals of Ioninfo structure to the log file
  dft_log("\n");
  print(dft_global_log, e);
  dft_log("\n");
}

// Print internal info for each species
void Ioninfo::print(Output *out, Everything &e)
{
  out->printf("Displaying Ioninfo contents:\n");
  out->printf("Found %d species:\n",nspecies);
  for (int sp = 0; sp < nspecies; sp++)
    species[sp].print(out, e);
}


/*
 * Make sure no atoms are layed on top of each other.
 */
void Ioninfo::check_pos()
{
  int dieflag = 0;
  real sizetest = 0;
  vector3 vtest[2];

  for (int sp=0; sp < nspecies; sp++)
    for (int n=0; n < species[sp].natoms; n++) {
      vtest[0] = species[sp].atpos[n];
      for (int sp1 = sp; sp1 < nspecies; sp1++)
	for (int n1 = ((sp1==sp) ? (n+1) : 0);
	     n1 < species[sp1].natoms; n1++) {
	  vtest[1] = vtest[0] - species[sp1].atpos[n1];
	  sizetest = vtest[1] * vtest[1];
	  if (sizetest < MIN_ION_DISTANCE) {
	    dieflag = 1;
	    dft_log(">>Atoms' positions coincide!!\n");
	    dft_log(">>(sp %d : at %d) = (sp %d : at %d) : ",
		       sp, n, sp1, n1);
	    species[sp].atpos[n].print(dft_global_log," %f ");
	  }
	}
    }

  if(dieflag != 0)
    die(">> You put atoms on top of each other, moron. I have to quit now. Bye\n");
}


int Ioninfo::symmetrize_force(const Symmetries &symm)
{
  if (symm.nrot <= 1) // no need to do any symmetrization
    return 0;

  // for each atom, each symmetry maps it to one of its own species.

  for (int sp = 0; sp < nspecies; sp++)
    {
      vector3 *forces_tmp;
      int iat;

      int natoms = species[sp].natoms;
      vector3 *forces = species[sp].forces;
      forces_tmp = (vector3*)mymalloc(sizeof(vector3)*natoms,
				      "forces_tmp",
				      "symmetrize_force");
      for (iat = 0; iat < natoms; iat++)
	forces_tmp[iat] = 0;
      for (iat = 0; iat < natoms; iat++)
	for (int irot = 0; irot < symm.nrot; irot++)
	  forces_tmp[symm.maps[sp][irot][iat]] +=
	                      symm.f_sym[irot] * forces[iat];
      for (iat = 0; iat < natoms; iat++)
	forces[iat] = forces_tmp[iat] / symm.nrot;
      myfree(forces_tmp);
    }
  return 1;

}

void Ioninfo::print_force(const Output *out, Everything &e) const
{
  for (int sp=0; sp < nspecies; sp++)
    species[sp].print_force(out, e);
}