#include <stdio.h>
#include <math.h>
#include <string.h>
#include <objects.h>

#ifdef DOUBLE
#include <dutil.h>
#include <dblock.h>
#define pi_over_180 (atan((double) 1.0) / (double) 45.0)
#define float double
#include <dfragment.h>
#include <dfield.h>
#endif
#
#ifndef DOUBLE
#include <util.h>
#include <block.h>
#define pi_over_180 ( (float ) (atan((double) 1.0) / (double) 45.0) )
#include <fragment.h>
#include <field.h>
#endif

#include "tcl.h"

#define debug_xtal 0

#include <chemsh.h>


#ifdef NOC_
#define make_cluster_  make_cluster
#define periodic_potential_  periodic_potential
#define cluster_potential_ cluster_potential
#define load_asymmetric_unit_ load_asymmetric_unit
#define store_cluster_ store_cluster
#define init_cluster_build_ init_cluster_build
#define end_cluster_build_ end_cluster_build
#define append_cluster_atom_ append_cluster_atom
#define fit_potential_ fit_potential
#define load_xtal_data_ load_xtal_data
#endif


#ifdef LINUXF2C
#define make_cluster_  make_cluster__
#define periodic_potential_  periodic_potential__
#define cluster_potential_ cluster_potential__
#define load_asymmetric_unit_ load_asymmetric_unit__
#define store_cluster_ store_cluster__
#define init_cluster_build_ init_cluster_build__
#define end_cluster_build_ end_cluster_build__
#define append_cluster_atom_ append_cluster_atom__
#define fit_potential_ fit_potential__
#define load_xtal_data_ load_xtal_data__
#endif


#ifdef CRAYXX
#define make_cluster_         MAKE_CLUSTER
#define periodic_potential_   PERIODIC_POTENTIAL
#define cluster_potential_    CLUSTER_POTENTIAL
#define load_asymmetric_unit_ LOAD_ASYMMETRIC_UNIT
#define store_cluster_        STORE_CLUSTER
#define init_cluster_build_   INIT_CLUSTER_BUILD
#define end_cluster_build_    END_CLUSTER_BUILD
#define append_cluster_atom_  APPEND_CLUSTER_ATOM
#define fit_potential_        FIT_POTENTIAL
#define load_xtal_data_       LOAD_XTAL_DATA
#endif


static Frag xtal_frag;

/*
FRAG_make_prim
Construct unit cell contents
*/

FRAG_make_prim(Frag f)
{
  INT ierr, mode, nseed, iseed[10], iprint, ndim;
  double rcut;
  xtal_frag = f;
  mode = 0;
  iprint=0;
  ndim = f->ndim;
  make_cluster_(&ndim, &mode, &nseed, iseed, &rcut, &iprint, &ierr);
  return (int) ierr;
}


FRAG_cut_cluster(Frag f, 
                 int iseed_in, 
                 double rcut)
{
  INT ierr, mode, nseed, iseed[10], iprint, ndim;
  xtal_frag = f;
  mode  = 1;
  nseed = 1;
  iprint=0;
  iseed[0] = iseed_in;
  ndim = f->ndim;
  make_cluster_(&ndim, &mode, &nseed, iseed, &rcut, &iprint, &ierr);

  /* 
     Now flag the result as non-periodic 
     Maybe in future it is better to keep the dimensionality
   */

  f->ndim = 0;

  return (int) ierr;
}


FRAG_find_molecule(Frag f, 
		   int *iseed,
		   int nseed)
{
  INT ierr, mode, iprint, ndim;
  double rcut;
  int i;

  xtal_frag = f;
  mode  = 2;
  iprint  = 0;
  ndim = f->ndim;

  printf("nseed %d ",nseed);
  for(i=0;i<nseed;i++)printf("%d ",iseed[i]);
  printf("\n");

  make_cluster_(&ndim, &mode, &nseed, iseed, &rcut, &iprint, &ierr);

  /* 
     Now flag the result as non-periodic 
     Maybe in future it is better to keep the dimensionality
   */

  f->ndim = 0;

  return (int) ierr;
}

/*

Compute periodic electrostatic potential on a grid.
Atomic charges must have been defined for the
atoms in the assymetric unit.

*/
FRAG_periodic_potential(Frag frag, 
			Field grid,
			Field data)
{
  INT ierr, ndim, ng;
  ObjList lgrid2;
  Field grid2;

  xtal_frag = frag;

  /* Convert field to irregular form if required */
  if(grid->regular ==0 || grid->regular == 1){

    lgrid2 = get_objlist("pote.temp","field", CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
    if(!lgrid2)return -1;
    grid2=(Field)(lgrid2->object->data);
    convert_to_irreg(grid, grid2);

  }else{
    grid2  = grid;
  }

  FIELD_copy(data, grid);
  FIELD_allocate_regular_memory(data, 1);
  ndim = frag->ndim;
  ng = grid2->ng,
  periodic_potential_(&ndim, grid2->grid, data->data, &ng, &ierr);

  /* store the fact that it is scalar data */
  data->no = 1;

  /* Define the unit */

  if( grid != grid2) rel_objlist(lgrid2);

  return (int) ierr;
}

/*doc
Compute cluster electrostatic potential on a grid.
Atomic charges must have been defined.
*/

FRAG_cluster_potential(Frag frag, 
		       Field grid,
		       Field data)
{
  INT ierr, ng, na;
  ObjList lgrid2;
  Field grid2;
  double *c, *chg;
  ierr=0;


  /* Convert field to irregular form if required */
  if(grid->regular ==0 || grid->regular == 1){

    lgrid2 = get_objlist("cluster_potential.temp","field", CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
    if(!lgrid2)return -1;
    grid2=(Field)(lgrid2->object->data);
    convert_to_irreg(grid, grid2);

  }else{
    grid2  = grid;
  }


  FIELD_copy(data, grid);
  FIELD_allocate_regular_memory(data, 1);

  na = (INT) FRAG_get_number_of_centres(frag);
  c =   (double *) ckalloc(sizeof(double) * 3 * na );
  chg = (double *) ckalloc(sizeof(double) * na );
  FRAG_get_abs_coordinates_and_charges(frag, c, chg);

  ng = grid2->ng,
  cluster_potential_(grid2->grid, data->data, &ng, c, chg, &na);

  ckfree( (char *) chg);
  ckfree( (char *) c);

  /* store the fact that it is scalar data*/
  data->no = 1;

  if( grid != grid2) rel_objlist(lgrid2);

  return (int) ierr;
}

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

 Load xtal data - called from the fortran code to load CRYSTAL data
 structures from fragment object.

************************************************************************/
void load_xtal_data_(INT *ngrmps, 
                     INT *iflags, 
                     STRING_ARGS1,
                     INT *itnsg,
                     double *a1,
		     double *alpha, 
		     double *beta, 
		     double *gamma,
		     INT *xnat, 
                     INT *xnbq, 
                     INT *xnsh,
		     INT *ierr)
{
  int i, iret;
  char *hmssg;
  Frag f;
  int l;
  double cnv_len;

  f = xtal_frag;

  *ierr=0;

  hmssg = string1;

  if(!strcmp(f->space_group_sym,"????????????????") && itnsg == 0){
    printf("No space group data\n");
    *ierr=-1;
    return;
  }

  /* Load up input symbol */
  l=strlen(f->space_group_sym);
  for(i=0;i<l;i++)hmssg[i]=f->space_group_sym[i];
  for(i=l;i<16;i++)hmssg[i]=' ';

  *itnsg = f->space_group;

  iret = conversion_factor("au","angstrom",&cnv_len);  
  if(iret){  *ierr=-1;  return; }

  *alpha=90.0;
  *beta=90.0;
  *gamma=90.0;
  a1[0]=1.0;
  a1[1]=1.0;
  a1[2]=1.0;

  switch (f->ndim){

  case 2:
    /* Load cell constants */
    a1[1]  = f->b * cnv_len;
    *gamma = f->gamma;
  case 1:
    a1[0] = f->a  * cnv_len;
  case 0:
    *ngrmps = f->space_group;
    break;

  case 3:
    /* Load cell constants */
    a1[0] = f->a  * cnv_len;
    a1[1] = f->b  * cnv_len;
    a1[2] = f->c  * cnv_len;
    *alpha = f->alpha;
    *beta = f->beta;
    *gamma = f->gamma;
    break;

  }
  *xnat = (INT) f->nxatom;
  *xnbq = (INT) f->nxbq;
  *xnsh = (INT) f->nxshell;
}

void
load_asymmetric_unit_(INT *iatn,
                      double *aicfmu,
                      double *ucha0,
                      INT *ierr){
  int ii,i;
  Frag f;
  float cnv, zcnv, ycnv;

  f = xtal_frag;
  *ierr=0;


  /* 
     Internal units are fractional with some angstroms
     for 1 and 2d systems.

     Here is the internal unit used to get back to
     aus in the crystal code 
  */

  cnv = 0.529177;
  ycnv = zcnv = 1.0;

  if(f->ndim < 3)zcnv=cnv;
  if(f->ndim < 2)ycnv=cnv;

  i=0;
  for(ii=0;ii<f->nxatom;ii++){
    iatn[i] = f->xatoms[ii].znum;
    aicfmu[3*i+0]=f->xatoms[ii].pos.x[0];
    aicfmu[3*i+1]=f->xatoms[ii].pos.x[1]*ycnv;
    aicfmu[3*i+2]=f->xatoms[ii].pos.x[2]*zcnv;
    ucha0[i] = f->xatoms[ii].charge;
    i++;
  }

  for(ii=0;ii<f->nxbq;ii++){
    iatn[i] = 0;
    aicfmu[3*i+0]=f->xbqs[ii].pos.x[0];
    aicfmu[3*i+1]=f->xbqs[ii].pos.x[1]*ycnv;
    aicfmu[3*i+2]=f->xbqs[ii].pos.x[2]*zcnv;
    ucha0[i] = f->xbqs[ii].charge;
    i++;
  }

  for(ii=0;ii<f->nxshell;ii++){
    iatn[i] = 0;
    aicfmu[3*i+0]=f->xshells[ii].pos.x[0];
    aicfmu[3*i+1]=f->xshells[ii].pos.x[1]*ycnv;
    aicfmu[3*i+2]=f->xshells[ii].pos.x[2]*zcnv;
    ucha0[i] = f->xshells[ii].charge;
    i++;
  }
}
                     

/*

   imap identifies the atom in the fractional coordinate system
   that the atom is generated from

*/
void store_cluster_(INT *pnat, 
		 INT *pnbq, 
		 INT *pnsh,
		 double *xa,
		 double *ucha,
		 INT *ndim,
		 double paret[3][3],
		 INT *imap,
		 STRING_ARGS1,
		 INT *ierr)
{
  int i, imcon;
  Frag f;
  double x, y, z, ch;
  double cell[9];
  int iz, map;
  char sym[20];
  char *hmssg;
  int parent;

  f = xtal_frag;

  *ierr=-1;

  if(debug_xtal){
    printf("Map arrays At: ");
    for(i=0;i<*pnat;i++){
      printf("%d ",imap[i]);
    }
    printf("BQ: ");
    for(i=0;i<*pnbq;i++){
      printf("%d ",imap[i+*pnat]);
    }
    printf("SH: ");
    for(i=0;i<*pnsh;i++){
      printf("%d ",imap[i+*pnat+*pnbq]);
    }
    printf("\n");
  }

  FRAG_allocate_atoms(f,*pnat);
  for(i=0;i<*pnat;i++){
    map=imap[i];
    if(FRAG_get_xatom(f, map, sym, &x, &y, &z, &iz))return;
    x = *xa++;
    y = *xa++;
    z = *xa++;
    ch = *ucha++;
    if(FRAG_set_atom(f, 1+i, sym, x, y, z))return;
    if(FRAG_set_atom_map(f, 1+i, map))return;
    if(f->atch)if(FRAG_set_atom_charge(f, 1+i, ch))return;
  }

  FRAG_allocate_bqs(f,(int) *pnbq);

  for(i=0;i<*pnbq;i++){
    map=imap[i + *pnat];
    if(FRAG_get_xbq(f, map, sym, &x, &y, &z, &ch))return;
    x = *xa++;
    y = *xa++;
    z = *xa++;
    ch = *ucha++;
    map=imap[i+*pnat];
    if(FRAG_set_bq(f, 1+i, sym, x, y, z, ch))return;
    if(FRAG_set_bq_map(f, 1+i, map))return;
  }

  FRAG_allocate_shells(f,(int) *pnsh);
  
  for(i=0;i<*pnsh;i++){
    map=imap[i+*pnat+*pnbq];
    if(FRAG_get_xshell(f, map, sym, &x, &y, &z, &ch, &parent))return;
    x = *xa++;
    y = *xa++;
    z = *xa++;
    ch = *ucha++;
    /* Parent will be either not present or incorrect, so must
       store directly here */
    if(FRAG_set_shell(f, 1+i, sym, x, y, z, ch, parent))return;
    if(FRAG_set_shell_map(f, 1+i, map))return;
  }

  /* Store space group symbol in case it has been updated */

  hmssg = string1;

  for(i=0;i<16;i++)f->space_group_sym[i]=hmssg[i];
  f->space_group_sym[16]='\0';

  /* DL_POLY cell is transpose of crystal paret */

  cell[0] = paret[0][0];
  cell[1] = paret[1][0];
  cell[2] = paret[2][0];
  cell[3] = paret[0][1];
  cell[4] = paret[1][1];
  cell[5] = paret[2][1];
  cell[6] = paret[0][2];
  cell[7] = paret[1][2];
  cell[8] = paret[2][2];

  switch(f->ndim){
  case 0:
    imcon=0;
    break;
  case 1:
    printf("Bad periodicity, cannot store a 1D system\n");
    return;
  case 2:
    imcon=6;
    break;
  case 3:
    imcon=3;
    break;
  }
  FRAG_set_cell(f, cell, imcon);

  *ierr=0;
}


struct tmp_atom_struct {
  char *symbol;
  double x,y,z;
  int map;
  int type;
  struct tmp_atom_struct *next, *tail;
};

static struct tmp_atom_struct *tmp_atoms;

void init_cluster_build_()
{
  tmp_atoms = NULL;
}

#define ATOM  1
#define BQ    2
#define SHELL 3

void end_cluster_build_(INT *ierr)
{
  char sym[16];
  int nat, nbq, nsh, index, index2, iz;
  struct tmp_atom_struct *a, *atmp;
  double x, y, z, chg;
  Frag f;
  int parent;

  *ierr=1;

  f = xtal_frag;

  /* 
     First establish the numbers of atoms/shells/bqs in 
     the cluster
  */

  nat=nbq=nsh=0;

  if(debug_xtal)printf(" counts %d %d %d\n",  f->nxatom,  f->nxbq, f->nxshell);

  for(a=tmp_atoms;a;a=atmp){
    atmp = a->next;

    if(debug_xtal)printf(" map %d\n", a->map );

    a->type = ATOM;

    if(a->map > f->nxatom){
      a->map -= f->nxatom; 
      a->type = BQ;

      if(a->map > f->nxbq){
	a->map -= f->nxbq; 
	a->type = SHELL;
      }
    }

    switch(a->type){
    case ATOM:  nat++; break;
    case BQ:    nbq++; break;
    case SHELL: nsh++; break;
    }
  }

  if(debug_xtal)printf(" counts %d %d %d\n",  nat, nbq, nsh);

  if(FRAG_allocate_atoms(f,nat))return;
  if(FRAG_allocate_bqs(f,nbq))return;
  if(FRAG_allocate_shells(f,nsh))return;

  index=1;
  for(a=tmp_atoms;a;a=atmp){
    atmp = a->next;

    switch(a->type){

    case ATOM:
      index2 = index;
      if(FRAG_get_xatom(f, a->map, sym, &x, &y, &z, &iz))return;
      if(FRAG_set_atom(f,index,sym,a->x,a->y,a->z))return;
      if(FRAG_set_atom_map(f,index2,a->map))return;
      if(f->atch){
	if(FRAG_get_xatom_charge(f,a->map,&chg))return;
	if(FRAG_set_atom_charge(f,index2,chg))return;
      }
      index++;
      break;

    case BQ:
      index2 = index - f->natom;
      if(FRAG_get_xbq(f, a->map, sym, &x, &y, &z, &chg))return;
      if(FRAG_set_bq(f,index2,sym,a->x,a->y,a->z,chg))return;
      if(FRAG_set_bq_map(f,index2,a->map))return;
      index++;
      break;

    case SHELL:
      index2 = index - f->natom - f->nbq;
      if(FRAG_get_xshell(f, a->map, sym, &x, &y, &z, &chg, &parent))return;
      /* parent will be either not present or incorrect, so must 
         store directly here */
      if(FRAG_set_shell(f,index2,sym,a->x,a->y,a->z, chg, parent))return;
      if(FRAG_set_shell_map(f,index2,a->map))return;
      index++;
      break;
    }
    ckfree((char *) a);
  }
  *ierr = 0;
}

void append_cluster_atom_(INT *map,
			double *x,
			double *y,
			double *z)
{
  struct tmp_atom_struct *head, *tail, *new;
  new = (struct tmp_atom_struct *) ckalloc(sizeof(struct tmp_atom_struct));
  if(tmp_atoms){
    head = tmp_atoms;
    tail = head->tail;
    tail->next = new;
    head->tail = new;
  } else {
    tmp_atoms = new;
    new->tail = new;
  }
  new->next = NULL;
  new->map = *map;
  new->x = *x;
  new->y = *y;
  new->z = *z;
}

/*
   values of flag

   0 : simple fit of values on fitting grid to values 
       on sample grid, no periodicity

   1 : fit charges at the atomic sites

   2 : fit charges at the periodic crystallographic sites (using xtal symmetry)
        - note that in this case the cluster must be a subset of the crystal
	Note that if a molcule contains bqs, these are fitted

*/
int FRAG_potfit(Frag frag,
		Field fit,
		Field sample,
		int ozer,
		int ocharg,     /* Apply charge constraint */
		double charge,  /* Charge constraint */
		int ofitsy,     /* Use symmetry equivalences */
		int oper,       /* Treat fitting array as periodic */
		int debug, 
		int verbose)
{

  int oatoms;
  ObjList lfit2, lsample2;
  Field fit2, sample2;
  int ncent, i, j;
  int *tmp;

  double *result;
  double *fitg;

  INT *symmetry;
  INT ierr;
  INT ndim, ng, ng2;
  INT *isymtab, *neq;
  INT maxsym, maxunq;
  INT *mask;
  INT ozer1, oper1, ocharg1 ,ofitsy1;


  /* We may need this to get space group symmetry if ocell && osym,  or oper */

  if(frag){
    xtal_frag = frag;
  }

  oatoms = (!fit);

  /* Convert fields to irregular form if required */
  if( oatoms ){

    /* Computing charges at atom sites 

       It is assumed that charges should be placed at all atom/shell/bq positions 

     Symmetry equivalence 

       It is assumed that atoms that originate from the same
       atom in the asymmetric unit will be treated as 
       symmetry related .. i.e. the symmetry of the atom 
       charge array is not used

       This way generating charges on the asymmetric unit is
       always unambigous

     */

    maxunq = ng = ncent = frag->natom + frag->nbq + frag->nshell;

    symmetry = (INT *) ckalloc(ng * sizeof(INT)) ;
    result = (double *) ckalloc(sizeof(double) * ncent);
    fitg = (double *) ckalloc(sizeof(double) * 3 * ncent);

    FRAG_get_abs_coordinate_matrix(frag, fitg);

    if(ofitsy && frag->space_group){

      /* Build symmetry relationships between the fitting centres
	 It is assumed that these are only present if the structure
	 has been generated from a crystal.

	 Hold the first real atom corresponding to each
	 fractional atom in tmp 
       */

      if(frag->nxatom){
	tmp = (int *) ckalloc(frag->nxatom * sizeof(int)) ;
	for(i=0;i<frag->nxatom;i++)tmp[i]=-1;
	for(i=0;i<frag->natom;i++){
	  if(tmp[frag->atoms[i].map] != -1){
	    symmetry[i] = tmp[frag->atoms[i].map]+1;
	  } else {
	    tmp[frag->atoms[i].map] = i;
	    symmetry[i] = i+1;
	  }
	}
	ckfree((char *) tmp);
      }

      if(frag->nxbq){
	tmp = (int *) ckalloc(frag->nxbq * sizeof(int)) ;
	for(i=0;i<frag->nxbq;i++)tmp[i]=-1;
	for(i=0;i<frag->nbq;i++){
	  if(tmp[frag->bqs[i].map] != -1){
	    symmetry[i + frag->natom] = tmp[frag->bqs[i].map]+1;
	  } else {
	    tmp[frag->bqs[i].map] = i;
	    symmetry[i + frag->natom] = i+1;
	  }
	}
	ckfree((char *) tmp);
      }

      if(frag->nxshell){
	tmp = (int *) ckalloc(frag->nxshell * sizeof(int)) ;
	for(i=0;i<frag->nxshell;i++)tmp[i]=-1;
	for(i=0;i<frag->nshell;i++){
	  if(tmp[frag->shells[i].map] != -1){
	    symmetry[i + frag->natom +frag->nbq ] = tmp[frag->shells[i].map]+1;
	  } else {
	    tmp[frag->shells[i].map] = i;
	    symmetry[i + frag->natom +frag->nbq ] = i+1;
	  }
	}
	ckfree((char *) tmp);
      }

      maxsym = 16;

    }else{
      for(i=0;i<ncent;i++)symmetry[i] = i+1;
      maxsym = 1;
    }

  } else {

    if(fit->regular ==0 || fit->regular == 1){
      lfit2 = get_objlist("pote.temp","field", CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
      if(!lfit2)return -1;
      fit2=(Field)(lfit2->object->data);
      convert_to_irreg(fit, fit2);
    }else{
      fit2  = fit;
    }
    /* We will return the charges as values on the grid */
    if(FIELD_allocate_regular_memory(fit,1))return -1;
    result = fit->data;
    fitg = fit2->grid;

    if(debug)printf("fit->symmetry %ld ofitsy %d\n",(long) fit->symmetry, ofitsy  );

    symmetry = (INT *) ckalloc(fit->ng * sizeof(INT)) ;

    if(fit->symmetry && ofitsy){

      for(i=0;i<fit2->ng;i++){
	symmetry[i] = (INT) fit->symmetry[i];
      }

      /* may need more intelligent estimate */
      maxsym = 16;
      maxunq = fit2->ng;

    } else {
      maxsym = 1;
      for(i=0;i<fit2->ng;i++)symmetry[i] = i+1;
      maxunq = fit2->ng;
    }

    ng = fit2->ng;

  }

  if(sample->regular ==0 || sample->regular == 1){
    lsample2 = get_objlist("pote.irreg.tmp","field", CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
    if(!lsample2)return -1;
    sample2=(Field)(lsample2->object->data);
    convert_to_irreg(sample, sample2);
  }else{
    sample2  = sample;
  }

  /* Create mask array */
  mask = (INT *) ckalloc(sample->ng * sizeof(INT)) ;

  if(sample->mask){
    for(i=0;i<sample->ng;i++){
      mask[i] = (INT) sample->mask[i];
    }
  }else{
    for(i=0;i<sample->ng;i++){
      mask[i] = 0;
    }
  }

  /* Create workspace for symmetry data */
  isymtab = (INT *) ckalloc(maxunq * maxsym * sizeof(INT));
  neq = (INT *) ckalloc(maxunq * sizeof(INT));

  /* Debug print */
  if(debug){
    printf("oatoms = %d ofitsy %d oper %d   maxunq %d maxsym %d \n",
	   oatoms,ofitsy,oper,maxunq, maxsym);
  }

  if(verbose){
    printf("Fitting grid\n");
    for(i=0;i<ng;i++){
      printf("%12.4f %12.4f %12.4f  %d\n", 
	     *(fitg+(3*i)),
	     *(fitg+(3*i+1)),
	     *(fitg+(3*i+2)),
	     symmetry[i] );
    }

    printf("Sampling grid\n");
    for(i=0;i<sample2->ng;i++){
      printf("%12.4f %12.4f %12.4f  %18.6f\n", 
	     *(sample2->grid+(3*i)),
	     *(sample2->grid+(3*i+1)),
	     *(sample2->grid+(3*i+2)),
	     *(sample2->data+i) );
    }
  }
  fflush(stdout);

  if(frag){
    ndim= (INT) frag->ndim;
  }else{
    ndim=0;
  }
  ofitsy1 = (INT) ofitsy;
  ng2 = (INT) sample2->ng;
  oper1 = (INT) oper;
  ozer1 = (INT) ozer;
  ocharg1 = (INT) ocharg;

  fit_potential_(&ndim, 
		 fitg, &ng, 
		 symmetry, &ofitsy1,  /* Fitting grid (with optional symmetry info) */
		 isymtab, &maxsym, &maxunq, neq,  /* Symmetry tables  (workspace) */
		  /* maxsym is the maximum number of equivalents*/
		  /* only needed if ofitsy = .true. dimension*/
		  /* (max # of eqivs) * (max # of uniq atoms) */
		 sample2->grid, sample->data, &ng2,    /* Sample data to fit to */
		 mask,                 /* mask array for the sampling grid */
		 result,               /* result charges */
		 &ocharg1, &charge,     /* optional charge contraint  */
		 &oper1,                /* treat the fitting grid as periodic */
		 &ozer1,                /* allow zero of potential to float */
		 &ierr );
  if(ierr){

    printf("Problem with PDC fit .. no charges stored\n");

  } else {

    /* Save the potential derived charges as atomic charges */
    if( oatoms ) {
      frag->atch = 1;
      j=0;
      for(i=0; i < frag->natom; i++){
	frag->atoms[i].charge = result[j++];
	frag->atoms[i].set_charge = 1;
      }
      for(i=0; i < frag->nbq;   i++)frag->bqs[i].charge   =result[j++];
      for(i=0; i < frag->nshell;i++)frag->shells[i].charge=result[j++];

      /* Set the charges on the asymmetric unit */
      if(frag->space_group && ofitsy)FRAG_set_asymmetric_charges(frag);
      ckfree((char *) result);
    }

  }

  if( !oatoms && (fit != fit2)) rel_objlist(lfit2);
  if( sample != sample2) rel_objlist(lsample2);

  ckfree((char *) mask);
  ckfree((char *) isymtab);
  ckfree((char *) neq);

  ckfree((char *) symmetry);

  return ierr;

}


int FRAG_delete_symmetry(Frag frag, 
			 Frag res)
{
  int t1[3];
  int t2[3];
  double c[9];
  double det;
  double x, y, z, xf, yf, zf;
  int i;

  FRAG_copy(res,frag);

  /* re-Allocate space for fractional coordinates */
  free_atom_list(res->xatoms, res->nxatom);
  res->nxatom = res->natom;
  res->xatoms = allocate_atom(res->nxatom);

  /* construct cartesian to fractional transformation matrix */

  if(res->ndim && res->space_group){

    switch (res->ndim){

    case 3:

      c[0+3*0] = res->pvec[0].x[0];
      c[0+3*1] = res->pvec[0].x[1];
      c[0+3*2] = res->pvec[0].x[2];

      c[1+3*0] = res->pvec[1].x[0];
      c[1+3*1] = res->pvec[1].x[1];
      c[1+3*2] = res->pvec[1].x[2];

      c[2+3*0] = res->pvec[2].x[0];
      c[2+3*1] = res->pvec[2].x[1];
      c[2+3*2] = res->pvec[2].x[2];

      det = invert_matrix(c,3,t1,t2);

      for(i=0;i<res->natom;i++){

	copy_frag_atom(res->xatoms,res->atoms,i,i);

	x = res->atoms[i].pos.x[0];
	y = res->atoms[i].pos.x[1];
	z = res->atoms[i].pos.x[2];

	xf = c[0+3*0]*x + c[1+3*0]*y + c[2+3*0]*z;
	yf = c[0+3*1]*x + c[1+3*1]*y + c[2+3*1]*z;
	zf = c[0+3*2]*x + c[1+3*2]*y + c[2+3*2]*z;

	res->xatoms[i].pos.x[0] = xf;
	res->xatoms[i].pos.x[1] = yf;
	res->xatoms[i].pos.x[2] = zf;

	if(FRAG_set_atom_map(res, 1+i, 1+i))return -1;

      }

      break;

    case 2:

      c[0+2*0] = res->pvec[0].x[0];
      c[0+2*1] = res->pvec[0].x[1];

      c[1+2*0] = res->pvec[1].x[0];
      c[1+2*1] = res->pvec[1].x[1];

      det = invert_matrix(c,2,t1,t2);

      for(i=0;i<res->natom;i++){

	copy_frag_atom(res->xatoms,res->atoms,i,i);

	x = res->atoms[i].pos.x[0];
	y = res->atoms[i].pos.x[1];
	z = res->atoms[i].pos.x[2];

	xf = c[0+2*0]*x + c[1+2*0]*y;
	yf = c[0+2*1]*x + c[1+2*1]*y;
	zf = z;

	res->xatoms[i].pos.x[0] = xf;
	res->xatoms[i].pos.x[1] = yf;
	res->xatoms[i].pos.x[2] = zf;

	if(FRAG_set_atom_map(res, 1+i, 1+i))return -1;
      }
    }
    res->space_group = 1;
    strcpy(res->space_group_sym,"p 1             ");
  }
  return 0;
}