#include "mmclib.h"
#include "getvalue.h"
#include "parse.h"
#include "lstsqr.h"

typedef int *pint;
typedef int **ppint;
typedef Real *pReal;
typedef Real **ppReal;
typedef Real ***pppReal;
typedef Complex *pComplex;

MultiMonteCarlo::MultiMonteCarlo(const Structure &_lattice, const Array<Array<int> > &_site_type_list, const iVector3d &_supercell, 
		const SpaceGroup &space_group, const LinkedList<MultiCluster> &cluster_list, 
		const Array<Array<Array<Real> > > &_corrfunc) :
               lattice(_lattice), site_type_list(_site_type_list), supercell(_supercell), corrfunc(_corrfunc), cur_rho(), cur_conc(), mu(), allowed_flip_site(),allowed_flip_before(),allowed_flip_after(), flip_span(-1,-1,-1) {
    Real max_clus_len=0.;
    {
      LinkedListIterator<MultiCluster> c(cluster_list);
      for ( ; c; c++) {
        Real l=get_cluster_length(c->clus);
        if (l>max_clus_len) max_clus_len=l;
      }
    }
    margin=find_sphere_bounding_box(lattice.cell,max_clus_len);
    iVector3d supercell_orig=supercell;
    for (int mult=1; 1 ; mult++) {
      supercell=mult*supercell_orig;
      int i;
      for (i=0; i<3; i++) {
        if (supercell(i) < margin(i)) break;
      }
      if (i==3) break;
    }
    total_box=supercell+2*margin;
    int spin_total_size=total_box(0)*total_box(1)*total_box(2)*lattice.atom_pos.get_size();
    spin=new SPIN_TYPE[spin_total_size];
    for (int i=0; i<spin_total_size; i++) {
      spin[i]=0;
    }
    cur_energy=0;
    int spin_size=supercell(0)*supercell(1)*supercell(2)*lattice.atom_pos.get_size();
    spin_orig=new SPIN_TYPE[spin_size];
    for (int i=0; i<spin_size; i++) {
      spin_orig[i]=0;
    }
    cur_disorder_param=0.;

    rspin_size=(Real)spin_size;

    site_in_cell=lattice.atom_pos.get_size();
    {
      Array<int> frozen_site(site_in_cell);
      zero_array(&frozen_site);
      active_site_in_cell=site_in_cell;
      if (file_exists("frozensite.in")) {
	ifstream file("frozensite.in");
	while (skip_delim(file)) {
	  int s=-1;
	  file >> s;
	  if (!frozen_site(s)) {active_site_in_cell--;}
	  frozen_site(s)=1;
	}
      }
      which_site=new int[active_site_in_cell];
      int s=0;
      for (int i=0; i<active_site_in_cell; i++,s++) {
	while (frozen_site(s)==1) {s++;}
	which_site[i]=s;
      }
    }
    nb_spin_val=new int[site_in_cell];
    for (int i=0; i<site_in_cell; i++) {
	nb_spin_val[i]=site_type_list(lattice.atom_type(i)).get_size();
    }
    nb_clusters=new int[site_in_cell];
    cluster_size=new pint[site_in_cell];
    site_offset=new ppint[site_in_cell];
    spin_val_clus=new pppReal[site_in_cell];
    eci=new pReal[site_in_cell];
    which_cluster=new pint[site_in_cell];
    total_clusters=cluster_list.get_size();
    cur_rho.resize(total_clusters);
    zero_array(&cur_rho);
    pcur_rho=cur_rho.get_buf();
    new_rho=new Real[total_clusters];
    rcluster_mult_per_atom=new Real[total_clusters];
    {
      LinkedListIterator<MultiCluster> c(cluster_list);
      for (int i=0; c; c++,i++) {
	rcluster_mult_per_atom[i]=calc_multiplicity(*c,lattice.cell,space_group.point_op,space_group.trans)/(Real)lattice.atom_pos.get_size();
      }
    }

    E_ref=0.;
    which_is_empty=-1;
    nb_point=0;
    {
      LinkedListIterator<MultiCluster> c(cluster_list);
      int cluster_index=0;
      for (; c; c++, cluster_index++) {
	  if (c->clus.get_size()==0) {which_is_empty=cluster_index;}
	  if (c->clus.get_size()==1) {nb_point++;}
	}
      which_is_point=new int[nb_point];
      point_mult=new Real[nb_point];
	c.init(cluster_list);
      cluster_index=0;
      int point_index=0;
      for (; c; c++, cluster_index++) {
            if (c->clus.get_size()==1) {
              which_is_point[point_index]=cluster_index;
              point_mult[point_index]=rcluster_mult_per_atom[cluster_index];
              point_index++;
            }
	}
    }
    cur_conc.resize(nb_point);
    mu.resize(nb_point);
    zero_array(&mu);
    pmu=mu.get_buf();

    rMatrix3d inv_cell=!lattice.cell;
    for (int s=0; s<lattice.atom_pos.get_size(); s++) {
      LinkedList<MultiCluster> long_cluster_list;
      LinkedList<int> long_cluster_i_list;
      LinkedListIterator<MultiCluster> c(cluster_list);
      int cluster_index=0;
      
      for ( ; c; c++,cluster_index++) {
	if (c->clus.get_size()!=0) {
            Array<MultiCluster> clusters;
            find_equivalent_clusters(&clusters, *c,space_group.cell,space_group.point_op,space_group.trans);
            for (int ec=0; ec<clusters.get_size(); ec++) {
	      for (int center_atom=0; center_atom<clusters(ec).clus.get_size(); center_atom++) {
		rVector3d lat_shift=lattice.atom_pos(s)-clusters(ec).clus(center_atom);
		if (is_int(inv_cell*lat_shift)) {
                  MultiCluster *to_add=new MultiCluster(clusters(ec).clus.get_size());
                  to_add->clus(0)=clusters(ec).clus(center_atom)+lat_shift;
		  to_add->site_type(0)=clusters(ec).site_type(center_atom);
		  to_add->func(0)=clusters(ec).func(center_atom);
                  int j=1;
                  for (int i=0; i<clusters(ec).clus.get_size(); i++) {
                    if (i!=center_atom) {
                      to_add->clus(j)=clusters(ec).clus(i)+lat_shift;
                      to_add->site_type(j)=clusters(ec).site_type(i);
                      to_add->func(j)=clusters(ec).func(i);
                      j++;
                    }
                  }
                  long_cluster_list << to_add;
                  long_cluster_i_list << new int(cluster_index);
                }
              }
            }
        }
      }
      nb_clusters[s]=long_cluster_list.get_size();
      cluster_size[s]=new int[nb_clusters[s]];
      eci[s]=new Real[nb_clusters[s]];
      which_cluster[s]=new int[nb_clusters[s]];
      site_offset[s]=new pint[nb_clusters[s]];
      spin_val_clus[s]=new ppReal[nb_clusters[s]];
      c.init(long_cluster_list);
      LinkedListIterator<int> c_i(long_cluster_i_list);
      for (int ic=0; ic<nb_clusters[s]; ic++, c++, c_i++) {
        cluster_size[s][ic]=c->clus.get_size();
	eci[s][ic]=0.;
        which_cluster[s][ic]=*c_i;
        site_offset[s][ic]=new int[c->clus.get_size()];
        spin_val_clus[s][ic]=new pReal[c->clus.get_size()];
        for (int j=0; j<c->clus.get_size(); j++) {
          int offset_in_cell=which_atom(lattice.atom_pos,c->clus(j),inv_cell);
          iVector3d offset_cell=to_int(inv_cell*(c->clus(j)-lattice.atom_pos(offset_in_cell)));
          site_offset[s][ic][j]=((offset_cell(0)*total_box(1) + offset_cell(1))*total_box(2) + offset_cell(2))*site_in_cell + offset_in_cell - s;
	  spin_val_clus[s][ic][j]=new Real[nb_spin_val[offset_in_cell]];
	  for (int l=0; l<nb_spin_val[offset_in_cell]; l++) {
            spin_val_clus[s][ic][j][l]=_corrfunc(c->site_type(j))(c->func(j))(l);
          }
        }
      }
    }
  }

MultiMonteCarlo::~MultiMonteCarlo(void) {
    for (int s=0; s<site_in_cell; s++) {
      for (int ic=0; ic<nb_clusters[s]; ic++) {
        for (int j=0; j<cluster_size[s][ic]; j++) {
	    delete[] spin_val_clus[s][ic][j];
        }
        delete[] spin_val_clus[s][ic];
        delete[] site_offset[s][ic];
      }
      delete[] spin_val_clus[s];
      delete[] site_offset[s];
      delete[] cluster_size[s];
      delete[] eci[s];
      delete[] which_cluster[s];
    }
    delete[] which_site;
    delete[] spin_val_clus;
    delete[] site_offset;
    delete[] cluster_size;
    delete[] eci;
    delete[] which_cluster;
    delete[] nb_clusters;

    delete[] nb_spin_val;
    delete[] new_rho;
    delete[] rcluster_mult_per_atom;
    delete[] which_is_point;
    delete[] point_mult;

    delete[] spin;
    delete[] spin_orig;
}

void MultiMonteCarlo::calc_delta_point_corr(Array<Real> *pcorr, int site, int type) {
  Array<Real> allcorr(total_clusters);
  zero_array(&allcorr);
  for (int c=0; c<nb_clusters[site]; c++) {
    if (cluster_size[site][c]==1) {
      allcorr(which_cluster[site][c])=spin_val_clus[site][c][0][type];
    }
  }
  for (int p=0; p<nb_point; p++) {
    (*pcorr)(p)+=allcorr(which_is_point[p]); // point_mult[p];
  }
}

void MultiMonteCarlo::calc_delta_point_corr(Array<Real> *pcorr, const Array<int> &sites,const Array<int> &types) {
  pcorr->resize(nb_point);
  zero_array(pcorr);
  for (int i=0; i<sites.get_size(); i++) {
    calc_delta_point_corr(pcorr,sites(i),types(i));
  }
}

void MultiMonteCarlo::find_all_allowed_flips(const Array2d<Real> &corr_constraints, const iVector3d &_flip_span, int min_nb_flip) {
  flip_span=_flip_span;
  int withinonecell=0;
  if (flip_span(0)==0 && flip_span(1)==0 && flip_span(2)==0) {withinonecell=1;}
  LinkedList<Array<int> > site_list;
  LinkedList<Array<int> > before_flip_list;
  LinkedList<Array<Array<int> > > after_flip_list;
  Array<int> site_done(site_in_cell);
  zero_array(&site_done);
  int nb_site_done=0;
  int nb_in_basis=0;
  Array<Array<Real> > corr_basis;
  int nbflip=1;
  while (1) {
    Array<int> maxsite(nbflip);
    for (int i=0; i<nbflip; i++) {maxsite(i)=active_site_in_cell;}
    MultiDimIterator<Array<int> > rawsites(maxsite);
    for ( ; rawsites; rawsites++) {
      Array<int> sites(maxsite.get_size());
      for (int i=0; i<sites.get_size(); i++) {
	sites(i)=which_site[((Array<int> &)rawsites)(i)];
      }
      Array<int> maxtype(maxsite.get_size());
      for (int i=0; i<maxsite.get_size(); i++) {maxtype(i)=nb_spin_val[((Array<int> &)sites)(i)];}
      MultiDimIterator<Array<int> > before_types(maxtype);
      int nb_tries=1;
      for (int i=0; i<maxtype.get_size(); i++) {nb_tries*=maxtype(i);}
      Array<int> before_flip(nb_tries);
      for (int i=0; i<before_flip.get_size(); i++) {before_flip(i)=-1;}
      int nb_flip=0;
      LinkedList<Array<int> > after_flip;
      for (; before_types; before_types++) {
        Array<Real> before_corr;
        calc_delta_point_corr(&before_corr, (const Array<int> &)sites,(const Array<int> &)before_types);
        MultiDimIterator<Array<int> > after_types(maxtype);
        for (; after_types; after_types++) {
          int same;
          for (same=0; same<nbflip; same++) {
            if (((Array<int> &)before_types)(same)==((Array<int> &)after_types)(same)) break;
	  }
	  if (same==nbflip) {
	    Array<Real> after_corr;
	    calc_delta_point_corr(&after_corr, (const Array<int> &)sites,(const Array<int> &)after_types);
	    Array<Real> delta_corr;
	    diff(&delta_corr,after_corr,before_corr);
	    Array<Real> c;
	    if (corr_constraints.get_size()(1)!=delta_corr.get_size()) {
	      ERRORQUIT("Incorrect number of point correlations to use the hybrid grandcanonical/canonical mode.\nPlease make sure clusters.out includes all point correlations");
	    }
	    product(&c,corr_constraints,delta_corr);
	    if (inner_product(c,c)<sqr(zero_tolerance)) {
	      int index=0;
	      for (int i=0; i<maxtype.get_size(); i++) {index=index*maxtype(i)+((Array<int> &)before_types)(i);}
	      before_flip(index)=nb_flip;
	      after_flip << new Array<int>(after_types);
	      nb_flip++;
	      if (min_nb_flip==0) {
		for (int l=0; l<((Array<int> &)sites).get_size();  l++) {
		  cerr << ((Array<int> &)sites)(l) << ": " << ((Array<int> &)before_types)(l) << "->" << ((Array<int> &)after_types)(l) << endl;
		}
		cerr << endl;
	      }
	      build_basis(&corr_basis,&nb_in_basis,delta_corr);
	    }
	  }
	}
      }
      if (nb_flip>0) {
	for (int i=0; i<((Array<int> &)sites).get_size(); i++) {
	  if (site_done(((Array<int> &)sites)(i))==0) {
	    site_done(((Array<int> &)sites)(i))=1;
	    nb_site_done++;
	  }
	}
	int addit=1;
	if (withinonecell) {
	  for (int i=0; i<sites.get_size(); i++) {
	      for (int j=0; j<i; j++) {
		if (sites(j)==sites(i)) {
		  addit=0;
		  break;
		}
	      }
	      if (!addit) {break;}
	  }
	}
	if (addit) {
	  site_list << new Array<int>(sites);
	  before_flip_list << new Array<int>(before_flip);
	  Array<Array<int> > *pafter_flip_array=new Array<Array<int> >();
	  LinkedList_to_Array(pafter_flip_array,after_flip);
	  after_flip_list << pafter_flip_array;
	}
	else {
	  cerr << "Ignored" << endl;
	}
      }
    }
    if (nb_in_basis >= nb_point-corr_constraints.get_size()(0) && (nb_site_done==active_site_in_cell) && nbflip>=min_nb_flip) break;
    nbflip++;
  }
  LinkedList_to_Array(&allowed_flip_site,site_list);
  LinkedList_to_Array(&allowed_flip_before,before_flip_list);
  LinkedList_to_Array(&allowed_flip_after,after_flip_list);
}

int MultiMonteCarlo::force_spin_flip(int *cell, int incell, int newspin) {
  int i;
  int mcell[3],mcellscan[3],offset;
  int oldspin;
  Real denergy,dconc;
  int cluster_count,site_count;
  int **pcluster;
  Real ***pppspin_val_clus,**ppspin_val_clus;
  int *pwhich_cluster;
  int *poffset,*psize;
  Real *peci, *prho, *prho2, *prcluster_mult_per_atom;
  Real rho;
  for (i=0; i<3; i++) {
    mcell[i]=cell[i]+margin(i);
  }

  offset=((mcell[0]*total_box(1) + mcell[1])*total_box(2) + mcell[2])*site_in_cell + incell;
  oldspin=spin[offset];
  denergy=0.;
  for (i=total_clusters, prho=new_rho; i>0; i--, prho++) {
    *prho=0.;
  }
  for (cluster_count=nb_clusters[incell], pcluster=site_offset[incell], pppspin_val_clus=spin_val_clus[incell], psize=cluster_size[incell], pwhich_cluster=which_cluster[incell], peci=eci[incell]; cluster_count>0; cluster_count--, pcluster++, pppspin_val_clus++, psize++, pwhich_cluster++, peci++) {
    poffset=*pcluster;
    ppspin_val_clus=*pppspin_val_clus;
    rho=(*ppspin_val_clus)[newspin]-(*ppspin_val_clus)[oldspin];
    poffset++;
    ppspin_val_clus++;
    for (site_count=(*psize)-1; site_count>0; site_count--, poffset++, ppspin_val_clus++) {
      rho*=(*ppspin_val_clus)[spin[offset+(*poffset)]];
    }
    new_rho[*pwhich_cluster]+=rho;
    denergy+=(*peci)*rho;
  }
  for (i=0; i<nb_point; i++) {
    denergy-=pmu[i]*new_rho[which_is_point[i]];
  }

  cur_energy+=denergy/rspin_size;
  for (i=0, prho=new_rho, prho2=pcur_rho, prcluster_mult_per_atom=rcluster_mult_per_atom; i<total_clusters; i++, prho++, prho2++, prcluster_mult_per_atom++) {
    *prho2+=(*prho)/(*prcluster_mult_per_atom)/rspin_size;
  }

  for (i=0; i<3; i++) {
    if (mcell[i]>=supercell(i)) {mcell[i]-=supercell(i);}
  }
  for (mcellscan[0]=mcell[0]; mcellscan[0]<total_box(0); mcellscan[0]+=supercell(0)) {
    for (mcellscan[1]=mcell[1]; mcellscan[1]<total_box(1); mcellscan[1]+=supercell(1)) {
      for (mcellscan[2]=mcell[2]; mcellscan[2]<total_box(2); mcellscan[2]+=supercell(2)) {
	spin[((mcellscan[0]*total_box(1) + mcellscan[1])*total_box(2) + mcellscan[2])*site_in_cell + incell]=newspin;
      }
    }
  }
  offset=((cell[0]*supercell(1) + cell[1])*supercell(2) + cell[2])*site_in_cell + incell;
  cur_disorder_param+=((newspin!=spin_orig[offset])-(oldspin!=spin_orig[offset]))/rspin_size;
  extension_update_spin_flip(cell,incell,oldspin,newspin);
  return oldspin;
}

void MultiMonteCarlo::save_state(MultiMonteCarloState *pstate) {
  pstate->cur_energy=cur_energy;
  pstate->cur_disorder_param=cur_disorder_param;
  pstate->cur_rho.resize(total_clusters);
  for (int i=0; i<total_clusters; i++) {pstate->cur_rho(i)=pcur_rho[i];}
  pstate->saved_spins.delete_all();
}

void MultiMonteCarlo::save_spin(MultiMonteCarloState *pstate, int *cell, int incell) {
  FlipInfo *pinfo=new FlipInfo();
  for (int i=0; i<3; i++) {
    pinfo->cell[i]=cell[i]+margin(i);
  }
  pinfo->incell=incell;
  int offset=((pinfo->cell[0]*total_box(1) + pinfo->cell[1])*total_box(2) + pinfo->cell[2])*site_in_cell + incell;
  pinfo->spin=spin[offset];
  pstate->saved_spins << pinfo;
}

void MultiMonteCarlo::restore_state(const MultiMonteCarloState &state) {
  cur_energy=state.cur_energy;
  cur_disorder_param=state.cur_disorder_param;
  for (int i=0; i<total_clusters; i++) {pcur_rho[i]=state.cur_rho(i);}
  LinkedListIterator<FlipInfo> it(state.saved_spins);
  for ( ; it; it++) {
    int mcell[3];
    for (int i=0; i<3; i++) {
      mcell[i]=it->cell[i];
      if (mcell[i]>=supercell(i)) {mcell[i]-=supercell(i);}
    }
    int mcellscan[3];
    for (mcellscan[0]=mcell[0]; mcellscan[0]<total_box(0); mcellscan[0]+=supercell(0)) {
      for (mcellscan[1]=mcell[1]; mcellscan[1]<total_box(1); mcellscan[1]+=supercell(1)) {
	for (mcellscan[2]=mcell[2]; mcellscan[2]<total_box(2); mcellscan[2]+=supercell(2)) {
	  spin[((mcellscan[0]*total_box(1) + mcellscan[1])*total_box(2) + mcellscan[2])*site_in_cell + it->incell]=it->spin;
	}
      }
    }
  }
}

int MultiMonteCarlo::access(int *cell, int incell) {
  int mcell[3];
  for (int i=0; i<3; i++) {
    mcell[i]=cell[i]+margin(i);
  }
  int offset=((mcell[0]*total_box(1) + mcell[1])*total_box(2) + mcell[2])*site_in_cell + incell;
  return spin[offset];
}

void MultiMonteCarlo::multi_spin_flip(void) {
  int which_flip,before_index,after_index;
  Array<iVector3d> cells;
  Array<int> oldspin;
  do {
    which_flip=random(allowed_flip_site.get_size());
    cells.resize(allowed_flip_site(which_flip).get_size());
    iVector3d center;
    if (flip_span(0)>=0) {
      for (int j=0; j<3; j++) {
	center(j)=random(supercell(j));
      }
    }
    for (int i=0; i<cells.get_size(); i++) {
      int ii;
      do {
	for (int j=0; j<3; j++) {
	  if (flip_span(0)>=0) {
	    cells(i)(j)=(supercell(j)+center(j)+random(2*flip_span(j)+1)-flip_span(j)) % supercell(j);
	  }
	  else {
	    cells(i)(j)=random(supercell(j));
	  }
	}
	for (ii=0; ii<i; ii++) {
	  if (allowed_flip_site(which_flip)(ii)==allowed_flip_site(which_flip)(i) && cells(ii)==cells(i)) break;
	}
      } while (ii!=i);
    }
    before_index=0;
    oldspin.resize(cells.get_size());
    for (int i=0; i<cells.get_size(); i++) {
      before_index=before_index*nb_spin_val[allowed_flip_site(which_flip)(i)]+access(cells(i).get_buf(),allowed_flip_site(which_flip)(i));
      oldspin(i)=access(cells(i).get_buf(),allowed_flip_site(which_flip)(i));
    }
    after_index=allowed_flip_before(which_flip)(before_index);
  } while (after_index==-1);
  MultiMonteCarloState saved_state;
  save_state(&saved_state);
  Real before_energy=cur_energy;
  Real before_extension_energy=extension_get_energy();
  extension_save_state();
  //  cerr << "---<" << endl;
  //  cerr << get_cur_energy() << " " << (-0.5*get_cur_concentration()(0)+1.0*get_cur_concentration()(1)) << endl;
  //  cerr << which_flip << endl << before_index << endl << after_index << endl;
  for (int i=0; i<cells.get_size(); i++) {
    save_spin(&saved_state,cells(i).get_buf(),allowed_flip_site(which_flip)(i));
    force_spin_flip(cells(i).get_buf(),allowed_flip_site(which_flip)(i),allowed_flip_after(which_flip)(after_index)(i));
    //extension_update_spin_flip(cells(i).get_buf(),allowed_flip_site(which_flip)(i),oldspin(i),allowed_flip_after(which_flip)(after_index)(i));

    // cerr << i << ") " << allowed_flip_site(which_flip)(i) << allowed_flip_after(which_flip)(after_index)(i) << endl;
  }
  Real denergy=(cur_energy-before_energy+extension_get_energy()-before_extension_energy)*rspin_size;
  // cerr << get_cur_energy() << " " << (-0.5*get_cur_concentration()(0)+1.0*get_cur_concentration()(1)) << endl;
  int accept=0;
  if (denergy<0) {
    accept=1;
  }
  else {
    if (uniform01()<exp(-denergy/T)) accept=1;
  }
  if (!accept) {
    restore_state(saved_state);
    extension_undo_spin_flip();
  }
  else {
    extension_forget_state();
  }
  // cerr << (accept ? 'A' : 'R') << endl;
  // cerr << get_cur_energy() << " " << (-0.5*get_cur_concentration()(0)+1.0*get_cur_concentration()(1)) << endl;
  // cerr << ">---" << endl;
}

void MultiMonteCarlo::set_eci(const Array<Real> &new_eci) {
  if (which_is_empty>=0) {
    E_ref=new_eci(which_is_empty)/lattice.atom_pos.get_size();
  } else {
    E_ref=0.;
  }
  for (int s=0; s<site_in_cell; s++) {
    for (int c=0; c<nb_clusters[s]; c++) {
      eci[s][c]=new_eci(which_cluster[s][c]);
    }
  }
  calc_from_scratch();
}

void MultiMonteCarlo::init_random(const Array<Array<Real> > &conc) {
    Array<Array<Real> > sconc;
    sconc=conc;
    for (int i=0; i<sconc.get_size(); i++) {
      Real s=0.;
      for (int j=0; j<sconc(i).get_size(); j++) {
         s+=conc(i)(j);
         sconc(i)(j)=s;
      }
      for (int j=0; j<sconc(i).get_size(); j++) {
         sconc(i)(j)/=s;
      }     
    }
    iMatrix3d msupercell;
    msupercell.diag(supercell);
    BoundingBox<int,3> bb(iVector3d(0,0,0),total_box-iVector3d(1,1,1));
    MultiDimIterator<iVector3d> cell(supercell);
    for ( ; cell; cell++) {
      for (int offset_in_cell=0; offset_in_cell<site_in_cell; offset_in_cell++) {
        Real r=uniform01();
        int the_spin=0;
        for (; the_spin<sconc(offset_in_cell).get_size()-1; the_spin++) {
          if (r<sconc(offset_in_cell)(the_spin)) break;
        }
        MultiDimIterator<iVector3d> image(iVector3d(-1,-1,-1),iVector3d(1,1,1));
        for ( ; image; image++) {
          iVector3d image_cell=(iVector3d &)cell+msupercell*(iVector3d &)image+margin;
          if (bb.is_in(image_cell)) {
            int offset=((image_cell(0)*total_box(1) + image_cell(1))*total_box(2) + image_cell(2))*site_in_cell + offset_in_cell;
            spin[offset]=the_spin;
          }
        }
      }
    }
    int spin_size=supercell(0)*supercell(1)*supercell(2)*site_in_cell;
    for (int i=0; i<spin_size; i++) {
      spin_orig[i]=0;
    }
    calc_from_scratch();
    extension_calc_from_scratch();
  }

void MultiMonteCarlo::init_structure(const Structure &str) {
  cerr << "begin init" << endl;
  int easycell=1;
  for (int i=0; i<3; i++) {
    int j;
    for (j=0; j<3; j++) {
      if (norm(lattice.cell.get_column(i)*(Real)supercell(i)-str.cell.get_column(j))<zero_tolerance) {break;}
    }
    if (j==3) {easycell=0;}
  }
  if (easycell) {
    cerr << "Easy cell detected!" << endl;
    iMatrix3d msupercell;
    msupercell.diag(supercell);
    BoundingBox<int,3> bb(iVector3d(0,0,0),total_box-iVector3d(1,1,1));
    for (int at=0; at<str.atom_pos.get_size(); at++) {
      rMatrix3d inv_cell=!(lattice.cell);
      int s=which_atom(lattice.atom_pos,str.atom_pos(at),inv_cell);
      iVector3d cell=to_int(inv_cell*(str.atom_pos(at)-lattice.atom_pos(s)));
      SPIN_TYPE newspin=(SPIN_TYPE)( str.atom_type(at)==-1 ? random(nb_spin_val[s]) : str.atom_type(at));
      MultiDimIterator<iVector3d> image(iVector3d(-1,-1,-1),iVector3d(1,1,1));
      for ( ; image; image++) {
	iVector3d image_cell=cell+msupercell*(iVector3d &)image+margin;
	if (bb.is_in(image_cell)) {
	  int offset=((image_cell(0)*total_box(1) + image_cell(1))*total_box(2) + image_cell(2))*site_in_cell + s;
	  spin[offset]=newspin;
	}
      }
      int offset=((cell(0)*supercell(1) + cell(1))*supercell(2) + cell(2))*site_in_cell + s;
      spin_orig[offset]=newspin;
    }
  }
  else {
    rMatrix3d inv_cell=!str.cell;
    iMatrix3d msupercell;
    msupercell.diag(supercell);
    BoundingBox<int,3> bb(iVector3d(0,0,0),total_box-iVector3d(1,1,1));
    MultiDimIterator<iVector3d> cell(supercell);
    int quickload=0;
    if (file_exists("quickload")) {quickload=1;}
    cerr << "quickload=" << quickload << endl;;
    int at=-1;
    for (; cell; cell++) {
      rVector3d rcell=lattice.cell*to_real((iVector3d &)cell);
      for (int s=0; s<site_in_cell; s++) {
	if (quickload) {
	  at++;
	}
	else {
	  at=which_atom(str.atom_pos,rcell+lattice.atom_pos(s),inv_cell);
	}
	SPIN_TYPE newspin=(SPIN_TYPE)( str.atom_type(at)==-1 ? random(nb_spin_val[s]) : str.atom_type(at));
	MultiDimIterator<iVector3d> image(iVector3d(-1,-1,-1),iVector3d(1,1,1));
	for ( ; image; image++) {
	  iVector3d image_cell=(iVector3d &)cell+msupercell*(iVector3d &)image+margin;
	  if (bb.is_in(image_cell)) {
	    int offset=((image_cell(0)*total_box(1) + image_cell(1))*total_box(2) + image_cell(2))*site_in_cell + s;
	    spin[offset]=newspin;
	  }
	}
	iVector3d thecell=(iVector3d &)cell;
	int offset=((thecell(0)*supercell(1) + thecell(1))*supercell(2) + thecell(2))*site_in_cell + s;
	spin_orig[offset]=newspin;
      }
    }
  }
  cerr << "end init" << endl;
  calc_from_scratch();
  extension_calc_from_scratch();
  cerr << "end calc" << endl;
}

void MultiMonteCarlo::calc_from_scratch(void) {
    cur_disorder_param=0.;
    cur_energy=0;
    for (int i=0; i<total_clusters; i++) {cur_rho(i)=0.;}
    MultiDimIterator<iVector3d> cur_cell(supercell);
    for ( ; cur_cell; cur_cell++) {
      iVector3d m_cur_cell=(iVector3d &)cur_cell+margin;
      for (int s=0; s<site_in_cell; s++) {
        int moffset=((m_cur_cell(0)*total_box(1) + m_cur_cell(1))*total_box(2) + m_cur_cell(2))*site_in_cell + s;
        iVector3d cell=cur_cell;
        int offset=((cell(0)*supercell(1) + cell(1))*supercell(2) + cell(2))*site_in_cell + s;
        cur_disorder_param+=(spin[moffset]!=spin_orig[offset]);
        for (int c=0; c<nb_clusters[s]; c++) {
          Real rho=1;
          for (int i=0; i<cluster_size[s][c]; i++) {
            rho*=spin_val_clus[s][c][i][spin[moffset+site_offset[s][c][i]]];
          }
	  rho/=(Real)(cluster_size[s][c]);
          cur_energy+=rho*eci[s][c];
	  pcur_rho[which_cluster[s][c]]+=rho/rcluster_mult_per_atom[which_cluster[s][c]];
        }
      }
    }
    cur_disorder_param/=rspin_size;
    cur_energy/=rspin_size;
    product(&cur_rho,cur_rho,1./rspin_size);
    if (which_is_empty>=0) {cur_rho(which_is_empty)=1.;}
    for (int i=0; i<nb_point; i++) {
      cur_energy-=mu(i)*cur_rho(which_is_point[i])*point_mult[i];
    }
  }

const Array<Real> & MultiMonteCarlo::get_cur_concentration(void) {
  for (int i=0; i<nb_point; i++) {
    cur_conc(i)=cur_rho(which_is_point[i])*point_mult[i];
  }
  return cur_conc;
}

void MultiMonteCarlo::set_T_mu(Real _T, const Array<Real> &_mu) {
    T=_T;
    for (int i=0; i<nb_point; i++) {
      cur_energy-=(_mu(i)-mu(i))*cur_rho(which_is_point[i])*point_mult[i];
    }
    mu=_mu;
    pmu=mu.get_buf();
}

void MultiMonteCarlo::spin_flip(void) {
    int i;
    int cell[3],mcell[3],mcellscan[3],incell,inactivecell,offset;
    int oldspin,newspin;
    Real denergy,dconc;
    int cluster_count,site_count;
    int **pcluster;
    Real ***pppspin_val_clus,**ppspin_val_clus;
    int *pwhich_cluster;
    int *poffset,*psize;
    Real *peci, *prho, *prho2, *prcluster_mult_per_atom;
    Real rho;
    int accept;
    for (i=0; i<3; i++) {
      cell[i]=random(supercell(i));
      mcell[i]=cell[i]+margin(i);
    }
    // incell=random(site_in_cell);
    inactivecell=random(active_site_in_cell);
    incell=which_site[inactivecell];
    offset=((mcell[0]*total_box(1) + mcell[1])*total_box(2) + mcell[2])*site_in_cell + incell;
    oldspin=spin[offset];
    newspin=(oldspin+1+random(nb_spin_val[incell]-1)) % nb_spin_val[incell];
    denergy=0.;
    for (i=total_clusters, prho=new_rho; i>0; i--, prho++) {
      *prho=0.;
    }
    for (cluster_count=nb_clusters[incell], pcluster=site_offset[incell], pppspin_val_clus=spin_val_clus[incell], psize=cluster_size[incell], pwhich_cluster=which_cluster[incell], peci=eci[incell]; cluster_count>0; cluster_count--, pcluster++, pppspin_val_clus++, psize++, pwhich_cluster++, peci++) {
      poffset=*pcluster;
      ppspin_val_clus=*pppspin_val_clus;
      rho=(*ppspin_val_clus)[newspin]-(*ppspin_val_clus)[oldspin];
      poffset++;
      ppspin_val_clus++;
      for (site_count=(*psize)-1; site_count>0; site_count--, poffset++, ppspin_val_clus++) {
        rho*=(*ppspin_val_clus)[spin[offset+(*poffset)]];
      }
      new_rho[*pwhich_cluster]+=rho;
      denergy+=(*peci)*rho;
    }
    for (i=0; i<nb_point; i++) {
      denergy-=pmu[i]*new_rho[which_is_point[i]];
    }
    extension_save_state();
    Real save_extension_energy=extension_get_energy();
    extension_update_spin_flip(cell,incell,oldspin,newspin);
    Real dextension_energy=(extension_get_energy()-save_extension_energy)*rspin_size;
    accept=0;
    if (denergy+dextension_energy<0) {
      accept=1;
    }
    else {
      if (uniform01()<exp(-(denergy+dextension_energy)/T)) accept=1;
    }
    if (accept) {
      cur_energy+=denergy/rspin_size;
      for (i=0, prho=new_rho, prho2=pcur_rho, prcluster_mult_per_atom=rcluster_mult_per_atom; i<total_clusters; i++, prho++, prho2++, prcluster_mult_per_atom++) {
         *prho2+=(*prho)/(*prcluster_mult_per_atom)/rspin_size;
      }

      for (i=0; i<3; i++) {
        if (mcell[i]>=supercell(i)) {mcell[i]-=supercell(i);}
      }
      for (mcellscan[0]=mcell[0]; mcellscan[0]<total_box(0); mcellscan[0]+=supercell(0)) {
        for (mcellscan[1]=mcell[1]; mcellscan[1]<total_box(1); mcellscan[1]+=supercell(1)) {
          for (mcellscan[2]=mcell[2]; mcellscan[2]<total_box(2); mcellscan[2]+=supercell(2)) {
            spin[((mcellscan[0]*total_box(1) + mcellscan[1])*total_box(2) + mcellscan[2])*site_in_cell + incell]=newspin;
          }
        }
      }
      offset=((cell[0]*supercell(1) + cell[1])*supercell(2) + cell[2])*site_in_cell + incell;
      cur_disorder_param+=((newspin!=spin_orig[offset])-(oldspin!=spin_orig[offset]))/rspin_size;
    }
    if (accept) {
      extension_forget_state();
    }
    else {
      extension_undo_spin_flip();
    }
}

void MultiMonteCarlo::run(int mc_passes, int mode) {
  int maxn=mc_passes*supercell(0)*supercell(1)*supercell(2)*site_in_cell;
  if (mode==1) {
    for (int n=0; n<maxn; n++) {
      spin_flip();
    }
  }
  else {
    for (int n=0; n<maxn; n++) {
      multi_spin_flip();
    }
  }
}

void MultiMonteCarlo::view(const Array<Arrayint> &labellookup, const Array<AutoString> &atom_label, ofstream &file, const rMatrix3d &axes) {
  for (int i=0; i<3; i++) {
    file << axes.get_column(i) << endl;
  }
  for (int i=0; i<3; i++) {
    file << (!axes)*((Real)(supercell(i))*lattice.cell.get_column(i)) << endl;
  }
  rMatrix3d iaxes=!axes;
  MultiDimIterator<iVector3d> cur_cell(supercell);
  for ( ; cur_cell; cur_cell++) {
    iVector3d m_cur_cell=(iVector3d &)cur_cell+margin;
    for (int s=0; s<site_in_cell; s++) {
      int moffset=((m_cur_cell(0)*total_box(1) + m_cur_cell(1))*total_box(2) + m_cur_cell(2))*site_in_cell + s;
      file << iaxes*(lattice.cell*to_real(cur_cell)+lattice.atom_pos(s)) << " " << atom_label(labellookup(lattice.atom_type(s))(spin[moffset])) << endl;
    }
  }
}

void MultiMonteCarlo::get_thermo_data(Array<Real> *pdata) {
  pdata->resize(2+nb_point+total_clusters);
  (*pdata)(0)=get_cur_energy();
  (*pdata)(1)=get_cur_disorder_param();
  int i=2;
  Array<Real> conc;
  conc=get_cur_concentration();
  for (int j=0; j<conc.get_size(); j++) {
    (*pdata)(i)=conc(j);
    i++;
  }
  Array<Real> rho;
  rho=get_cur_corr();
  for (int j=0; j<rho.get_size(); j++) {
    (*pdata)(i)=rho(j);
    i++;
  }
}

void run_mc(GenericAccumulator *accum, MultiMonteCarlo *pmc, int mode) {
  Array<Real> data;
  pmc->get_thermo_data(&data);
  while (accum->new_data(data)) {
    pmc->run(1,mode);
    pmc->get_thermo_data(&data);
  }
}

#include "anyfft.h"

KSpaceMultiMonteCarlo::KSpaceMultiMonteCarlo(const Structure &_lattice, const Array<Array<int> > &_site_type_list,
		const iVector3d &_supercell,
		const SpaceGroup &space_group, const LinkedList<MultiCluster> &cluster_list,
		const Array<Array<Array<Real> > > &_corrfunc, KSpaceECI *_p_kspace_eci):
  MultiMonteCarlo(_lattice,_site_type_list,_supercell,space_group,cluster_list,_corrfunc), flipped_spins(), ft_spin(), ft_eci(), dir_eci(), convol(), ref_x() {
  p_kspace_eci=_p_kspace_eci;
  nsite=_lattice.atom_pos.get_size();
  size=supercell(0)*supercell(1)*supercell(2);
  rsize=(Real)size;

  ft_spin.resize(nsite);
  for (int s=0; s<nsite; s++) {
    ft_spin(s).resize(nb_spin_val[s]);
    for (int ss=0; ss<ft_spin(s).get_size(); ss++) {
      ft_spin(s)(ss).resize(size);
    }
  }
  convol=ft_spin;

  ft_eci.resize(nsite);
  for (int s=0; s<nsite; s++) {
    ft_eci(s).resize(s+1);
    for (int t=0; t<=s; t++) {
      ft_eci(s)(t).resize(nb_spin_val[s]);
      for (int ss=0; ss<ft_eci(s)(t).get_size(); ss++) {
        ft_eci(s)(t)(ss).resize(nb_spin_val[t]);
        for (int tt=0; tt<ft_eci(s)(t)(ss).get_size(); tt++) {
          ft_eci(s)(t)(ss)(tt).resize(size);
        }
      }
    }
  }
  dir_eci=ft_eci;

  fft_time=0;
  flip_time=0;
  ref_x=get_cur_concentration();
  threshold_dx=1e-2;
  cur_recip_E=0.;
  psave_spins=NULL;
}

void plot_3d_surf(const rMatrix3d &k_mesh, const rMatrix3d &rec_lat, const iVector3d &supercell, Array<Complex> a) {
  {
    ofstream file("tmp.gnu");
    file << "set hidden" << endl 
	 << "splot 'tmp.out' u 1:2:3 w l" << endl
	 << "pause -1" << endl;
  }
    {
      ofstream file("tmp.out");
  for (int z=0; z<supercell(2); z++) {
      for (int y=0; y<supercell(1); y++) {
	for (int x=0; x<supercell(0); x++) {
	  int offset=((z*supercell(1) + y)*supercell(0) + x);
	      rVector3d k;
	      if (k_mesh==rec_lat) {
		 k=k_mesh*to_real(iVector3d(x,y,z));
	      } else {
		 k=flip_into_brillouin_1(k_mesh*to_real(iVector3d(x,y,z)),rec_lat);
	      }
	  file << k(0) << " " << k(1) << " " << k(2) << " " << real(a(offset)) << " " << imag(a(offset)) << endl;
	}
	file << endl;
      }
    }
    //    system("gnuplot tmp.gnu");
  }
}

void KSpaceMultiMonteCarlo::set_k_space_eci(void) {
  ref_x=get_cur_concentration();
  for (int s=0; s<nsite; s++) {
    for (int t=0; t<=s; t++) {
      for (int ss=0; ss<ft_eci(s)(t).get_size(); ss++) {
        for (int tt=0; tt<ft_eci(s)(t)(ss).get_size(); tt++) {
	  zero_array(&(ft_eci(s)(t)(ss)(tt)));
	}
      }
    }
  }
  p_kspace_eci->get_k_space_eci(&ft_eci, ref_x);
  dir_eci=ft_eci;
  for (int s=0; s<nsite; s++) {
    for (int t=0; t<=s; t++) {
      for (int ss=0; ss<dir_eci(s)(t).get_size(); ss++) {
        for (int tt=0; tt<dir_eci(s)(t)(ss).get_size(); tt++) {
          fftnd(dir_eci(s)(t)(ss)(tt).get_buf(),3,supercell.get_buf(),-1);
/*
rMatrix3d id;
id.identity();
plot_3d_surf(id,id,supercell,dir_eci(s)(t)(ss)(tt));
cerr << s << " " << t << " " << ss << " " << tt << " " << "eci" << endl;
*/
        }
      }
    }
  }
}

void KSpaceMultiMonteCarlo::extension_calc_from_scratch(void) {
//cerr << "b_total=" << E_ref+cur_energy+cur_recip_E+mu*cur_conc << endl;
  Array<Real> diffx;
  diff(&diffx,get_cur_concentration(),ref_x);
  if (max_norm(diffx)>threshold_dx) {
    set_k_space_eci();
  }
  fft_time=clock();
  flipped_spins.delete_all();
  psave_spins=NULL;
  flip_time=0;
  SPIN_TYPE *spin_corner=spin+((margin(0)*total_box(1) + margin(1))*total_box(2) + margin(2))*nsite;
  MultiDimIterator<iVector3d> cell(supercell);
  for ( ; cell; cell++) {
    iVector3d &rcell=(iVector3d &)cell;
    int offset=((rcell(0)*total_box(1) + rcell(1))*total_box(2) + rcell(2))*nsite;
    //int i=(rcell(0)*supercell(1) + rcell(1))*supercell(2) + rcell(2); //permbug
    int i=(rcell(2)*supercell(1) + rcell(1))*supercell(0) + rcell(0);
    for (int s=0; s<nsite; s++) {
      ft_spin(s)(0)(i)=Complex(1.,0.);
      for (int ss=1; ss<ft_spin(s).get_size(); ss++) {
        ft_spin(s)(ss)(i)=Complex(corrfunc(nb_spin_val[s]-2)(ss-1)(spin_corner[offset+s]),0.);
      }
    }
  }
//Array<Array<Array<Complex> > > savesp(ft_spin);
  for (int s=0; s<nsite; s++) {
    for (int ss=0; ss<ft_spin(s).get_size(); ss++) {
      fftnd(ft_spin(s)(ss).get_buf(),3,supercell.get_buf(),1);
//rMatrix3d id;
//id.identity();
//plot_3d_surf(id,id,supercell,ft_spin(s)(ss));
// cerr << "p" << endl;
    }
  }
  cur_recip_E=0.;
  for (int i=0; i<size; i++) {
    for (int s=0; s<nsite; s++) {
      for (int ss=0; ss<nb_spin_val[s]; ss++) {
        convol(s)(ss)(i)=0.;
        for (int t=0; t<nsite; t++) {
	  for (int tt=0; tt<nb_spin_val[t]; tt++) {
	      if (t<=s) {
		  convol(s)(ss)(i)+=ft_eci(s)(t)(ss)(tt)(i)*ft_spin(t)(tt)(i);
	      }
	      else {
		  convol(s)(ss)(i)+=conj(ft_eci(t)(s)(tt)(ss)(i))*ft_spin(t)(tt)(i);
	      }
          }
        }
      }
    }
  }
  for (int i=0; i<size; i++) {
    for (int s=0; s<nsite; s++) {
      for (int ss=0; ss<ft_spin(s).get_size(); ss++) {
	cur_recip_E+=real(conj(ft_spin(s)(ss)(i))*convol(s)(ss)(i));
//if (fabs(cur_recip_E)>0.) {
//cerr << cur_recip_E << endl;
//}
      }
    }
  }
//cerr << "end" << endl;
  for (int s=0; s<nsite; s++) {
    for (int ss=0; ss<convol(s).get_size(); ss++) {
      fftnd(convol(s)(ss).get_buf(),3,supercell.get_buf(),-1);
    }
  }
/*
  Real debE=0.;
  for (int s=0; s<nsite; s++) {
    for (int ss=0; ss<convol(s).get_size(); ss++) {
	for (int i=0; i<size; i++) {
	    debE+=real(convol(s)(ss)(i)*savesp(s)(ss)(i));
	}
    }
  }
  cerr << "debcE=" << debE/rspin_size << endl;

  {
      Real debE=0.;
      MultiDimIterator<iVector3d> cella(supercell);
      for ( ; cella; cella++) {
	  MultiDimIterator<iVector3d> cellb(supercell);
	  for ( ; cellb; cellb++) {
	      iVector3d &rcella=(iVector3d &)cella;
	      iVector3d &rcellb=(iVector3d &)cellb;
	      iVector3d dr=-(cellb-cella);
	      for (int j=0; j<3; j++) {
		  dr(j)=(dr(j)+supercell(j)) % supercell(j);
	      }
	      int offseta=((rcella(2)*supercell(1) + rcella(1))*supercell(0) + rcella(0));
	      int offsetb=((rcellb(2)*supercell(1) + rcellb(1))*supercell(0) + rcellb(0));
	      int offsetd=((dr(2)*supercell(1) + dr(1))*supercell(0) + dr(0));
	      int offsetp=((((supercell(2)-dr(2))%supercell(2))*supercell(1) + ((supercell(1)-dr(1))%supercell(1)))*supercell(0) + ((supercell(0)-dr(0))%supercell(0)));
	      for (int s=0; s<nsite; s++) {
		  for (int t=0; t<nsite; t++) {
		      for (int ss=0; ss<nb_spin_val[s]; ss++) {
			  for (int tt=0; tt<nb_spin_val[t]; tt++) {
			      if (t<=s) {
				  debE+=real(dir_eci(s)(t)(ss)(tt)(offsetd)*savesp(s)(ss)(offseta)*savesp(t)(tt)(offsetb));
			      } else {
				  debE+=real(dir_eci(t)(s)(tt)(ss)(offsetp)*savesp(s)(ss)(offseta)*savesp(t)(tt)(offsetb));
			      }
			  }
		      }
		  }
	      }
	  }
      }
      cerr << "debE=" << debE/rspin_size << endl;
  }
*/

  cur_recip_E/=(rsize*rspin_size);
cerr << "rE=" << cur_recip_E << endl;
  fft_time=clock()-fft_time;
}

void KSpaceMultiMonteCarlo::extension_save_state(void) {
  LinkedListIterator<FlippedSpin> tmp(flipped_spins);
  psave_spins=(FlippedSpin *)tmp;
  save_recip_E=cur_recip_E;
//cerr << "save" << flipped_spins.get_size() << endl;
}

void KSpaceMultiMonteCarlo::extension_forget_state(void) {
  psave_spins=NULL;
  int nbflip=flipped_spins.get_size();
  Array<Real> diffx;
  diff(&diffx,get_cur_concentration(),ref_x);
//cerr << "forget" << endl;
  if (max_norm(diffx)>threshold_dx) {
    extension_calc_from_scratch();
  } else {
    if (flip_time>0) {
//cerr << "time: " << nbflip << " " << flip_time << " " << 2*fft_time/flip_time << endl;
      if (nbflip > 2*fft_time/flip_time) {
//cerr << "before= " << extension_get_energy() << endl;
	  extension_calc_from_scratch();
//cerr << "after=  " << extension_get_energy() << endl;
      }
    }
  }
}

void KSpaceMultiMonteCarlo::extension_update_spin_flip(int *cell, int incell, int oldspin, int newspin) {

  flip_time=clock();
  int dr[3];
  int mdr[3];
  int *psupercell=supercell.get_buf();

  Array<Real> rspinchange(nb_spin_val[incell]);
  rspinchange(0)=0.;
  for (int ss=1; ss<rspinchange.get_size(); ss++) {
    rspinchange(ss)=corrfunc(rspinchange.get_size()-2)(ss-1)(newspin)-corrfunc(rspinchange.get_size()-2)(ss-1)(oldspin);
  }


  Real dE=0.;
  for (int ss=1; ss<nb_spin_val[incell]; ss++) {
    for (int tt=1; tt<nb_spin_val[incell]; tt++) {
      dE+=real(dir_eci(incell)(incell)(ss)(tt)(0)*rspinchange(ss)*rspinchange(tt));
    }
  }

//dE=0.;

  //int offset=((cell[0]*psupercell[1] + cell[1])*psupercell[2] + cell[2]); //permbug
  int offset=((cell[2]*psupercell[1] + cell[1])*psupercell[0] + cell[0]);

  for (int ss=1; ss<rspinchange.get_size(); ss++) {
    dE+=2.*real(convol(incell)(ss)(offset)*rspinchange(ss));
  }

  LinkedListIterator<FlippedSpin> i(flipped_spins);
  for (; i; i++) {
    for (int j=0; j<3; j++) {
      dr[j] =(cell[j]-(i->cell[j])+psupercell[j]) % psupercell[j];
      mdr[j]=(i->cell[j]-cell[j]+psupercell[j]) % psupercell[j];
    }
    //int offsetdr=((dr[0]*psupercell[1] + dr[1])*psupercell[2] + dr[2]); //permbug
    int offsetdr =((dr[2]*psupercell[1] + dr[1])*psupercell[0] + dr[0]);
    int offsetmdr=((mdr[2]*psupercell[1] + mdr[1])*psupercell[0] + mdr[0]);
    for (int ss=1; ss<nb_spin_val[incell]; ss++) {
      for (int tt=1; tt<nb_spin_val[i->incell]; tt++) {
        if (incell>=i->incell) {
          dE+=2.*real(dir_eci(incell)(i->incell)(ss)(tt)(offsetdr)*rspinchange(ss)*i->dspin(tt));
        } else {
          dE+=2.*real(dir_eci(i->incell)(incell)(tt)(ss)(offsetmdr)*rspinchange(ss)*i->dspin(tt));
	}
      }
    }
  }
  cur_recip_E+=dE/(rspin_size);

/*
 cerr << "Eafterflip=" << cur_recip_E << endl;
 int ck=0;
 if (ck==1) {
     extension_calc_from_scratch();
     cerr << cur_recip_E << endl;
 }
*/
  flipped_spins.push_front(new FlippedSpin(cell,incell,rspinchange));
  flip_time=clock()-flip_time;
}

void KSpaceMultiMonteCarlo::extension_undo_spin_flip(void) {
  LinkedListIterator<FlippedSpin> flip(flipped_spins);
  while ((FlippedSpin *)flip != psave_spins) {
    delete flipped_spins.detach(flip);
  }
  cur_recip_E=save_recip_E;
  psave_spins=NULL;
//cerr << "undo" << flipped_spins.get_size() << endl;
}

KSpaceMultiMonteCarlo::~KSpaceMultiMonteCarlo() {
  // free all dynamically allocated memory;
}