/*
    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 "command.h"

int ion_dep[] = {ION_SPECIES,LATTICE,COORDS_TYPE,0};

extern char processing_error[DFT_LINE_LEN];

int process_ion(char *input_line,Everything &e)
{
  char key[DFT_MSG_LEN];
  if (sscanf(input_line,"%*s %s",key)!=1)
    {
      sprintf(processing_error,
	      "ion command must be followed by a name.  Input line:\n%s\n",
	      input_line);
      return 1;
    }
  int sp;
  for (sp=0; sp < e.ioninfo.nspecies; sp++)
    if (MATCH(key,e.ioninfo.species[sp].name))
      {
	real x,y,z;
	int flag;
	if (sscanf(input_line,"%*s %*s %lg %lg %lg %d",&x,&y,&z,&flag)!=4)
	  {
	    sprintf(processing_error,
		    "ion must have x y z and move flag.  Input was:\n%s\n",
		    input_line);
	    return 1;
	  }
	if (flag != 0 && flag != 1)
	  {
	    sprintf(processing_error,
		    "move_flag for ion must be 1 or 0.  %d was input.\n",
		    flag);
	    return 1;
	  }
	int &nat = e.ioninfo.species[sp].natoms;
	nat++;
	Speciesinfo &spinfo = e.ioninfo.species[sp];
	if (nat==1)
	  {
	    spinfo.atpos = (vector3 *)mymalloc(sizeof(vector3)*nat,
					       "atpos[]","process_ion");
	    spinfo.forces = (vector3 *)mymalloc(sizeof(vector3)*nat,
						"forces[]","process_ion");
	    spinfo.old_forces = (vector3 *)mymalloc(sizeof(vector3)*nat,
						"old_forces[]","process_ion");
	    spinfo.ion_vel = (vector3 *)mymalloc(sizeof(vector3)*nat,
						   "ion_vel[]","process_ion");
	    spinfo.move = (int *)mymalloc(sizeof(int)*nat,
					  "move[]","process_ion");
	  }
	else
	  {
	    spinfo.atpos = (vector3 *)myrealloc(spinfo.atpos,
						sizeof(vector3)*nat,
						"atpos[]","process_ion");
	    spinfo.forces = (vector3 *)myrealloc(spinfo.forces,
						 sizeof(vector3)*nat,
						 "forces[]","process_ion");
	    spinfo.old_forces = (vector3 *)myrealloc(spinfo.old_forces,
						 sizeof(vector3)*nat,
						 "old_forces[]","process_ion");
	    spinfo.ion_vel = (vector3 *)myrealloc(spinfo.ion_vel,
						  sizeof(vector3)*nat,
						  "ion_vel[]","process_ion");
	    spinfo.move = (int *)myrealloc(spinfo.move,
					   sizeof(int)*nat,
					   "move[]","process_ion");
	  }
	if(e.ioninfo.coords_type == LATTICE_COORDS){
	  spinfo.atpos[nat-1].v[0] = x;
	  spinfo.atpos[nat-1].v[1] = y;
	  spinfo.atpos[nat-1].v[2] = z;
	} else if(e.ioninfo.coords_type == CARTESIAN_COORDS){
	  vector3 xx(x,y,z);

	  // HACK WARNING
	  //
	  // this is really stupid, but when the program gets here,
	  // the lattice is read in, but the lattice_scale has not been
	  // applied yet, so we need to apply it here.
	  //
	  // when this idiocy with the lattice_scale is sorted out,
	  // this might need to be changed!!!!!
	  matrix3 scaledR(e.lattice.latvec);
	  for (int i=0; i < 3; i++)
	    for (int j=0; j < 3; j++)
	      scaledR.m[j][i] *= e.cntrl.lattscale[i];


	  spinfo.atpos[nat-1] = inv3(scaledR)*xx;
	} else {
	  die("Invalid CoordsType in process_ion()\n");
	}

	spinfo.move[nat-1] = flag;
	return 0;
      }
  sprintf(processing_error,
	  "Species %s was not found.  Input line:\n%s\n",key,input_line);
  return 1;
}

void ion_print_status(Everything &e,Output *log)
{
  int sp,i;
  vector3 xx;
  for (sp=0; sp < e.ioninfo.nspecies; sp++)
    for (i=0; i < e.ioninfo.species[sp].natoms; i++){
      if(e.ioninfo.coords_type == LATTICE_COORDS){
	xx = e.ioninfo.species[sp].atpos[i];
      } else if(e.ioninfo.coords_type == CARTESIAN_COORDS){
	xx = e.lattice.latvec * e.ioninfo.species[sp].atpos[i];
      }
      log->printf("ion %s %18.14lf %18.14lf %18.14lf %d\n",
		  e.ioninfo.species[sp].name,
		  xx.v[0],
		  xx.v[1],
		  xx.v[2],
		  e.ioninfo.species[sp].move[i]);
    }
}

command* setup_ion(command *next)
{
  command* cmnd = alloc_command();

  cmnd->number = ION;
  cmnd->name = "ion";
  cmnd->format = "ion <id> <x1> <x2> <x3> <move_flag>";
  cmnd->comments = "# Specifies that ion type id at x y z.\n# move_flag is 1 | 0 and specifies whether the ion can be moved";
  cmnd->dependencies = ion_dep;
  cmnd->process_command = process_ion;
  cmnd->print_status = ion_print_status;
  cmnd->allow_multiple = 1;
  cmnd->next_command = next;

  return cmnd;
}