static int tclcommand_analyze_fluid_parse_densprof(Tcl_Interp *interp, int argc, char **argv) {
  
  int i, pdir, x1, x2;
  char buffer[TCL_DOUBLE_SPACE];
  double *profile;

  if (argc <3) {
    Tcl_AppendResult(interp, "usage: analyze fluid density <p_dir> <x1> <x2>", (char *)NULL);
    return TCL_ERROR;
  }

  if (!ARG_IS_I(0,pdir)) return TCL_ERROR;
  if (!ARG_IS_I(1,x1)) return TCL_ERROR;
  if (!ARG_IS_I(2,x2)) return TCL_ERROR;

  if (pdir != 2) {
    Tcl_AppendResult(interp, "analyze fluid density is only implemented for pdir=2 yet!", (char *)NULL);
    return TCL_ERROR;
  }

  profile = malloc(lblattice.grid[pdir]*node_grid[pdir]*sizeof(double));

  lb_master_calc_densprof(profile, pdir, x1, x2);

  for (i=0; i<lblattice.grid[pdir]*node_grid[pdir]; i++) {
    Tcl_PrintDouble(interp, i*lblattice.agrid, buffer);
    Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
    Tcl_PrintDouble(interp, profile[i], buffer);
    Tcl_AppendResult(interp, buffer, "\n", (char *)NULL);
  }
  free(profile);
  
  return TCL_OK;

}

/** Parser for the \ref lbnode_extforce command. Can be used in future to set more values like rho,u e.g.
*/
int tclcommand_lbnode_extforce_gpu(ClientData data, Tcl_Interp *interp, int argc, char **argv) {

  int err=TCL_ERROR;
  int coord[3];

  --argc; ++argv;
  
  if (argc < 3) {
    Tcl_AppendResult(interp, "too few arguments for lbnode_extforce", (char *)NULL);
    return TCL_ERROR;
  }

  if (!ARG_IS_I(0,coord[0]) || !ARG_IS_I(1,coord[1]) || !ARG_IS_I(2,coord[2])) {
    Tcl_AppendResult(interp, "wrong arguments for lbnode", (char *)NULL);
    return TCL_ERROR;
  } 
  argc-=3; argv+=3;

  if (argc == 0 ) { 
    Tcl_AppendResult(interp, "lbnode_extforce syntax: lbnode_extforce X Y Z [ print | set ] [ F(X) | F(Y) | F(Z) ]", (char *)NULL);
    return TCL_ERROR;
  }

  if (ARG0_IS_S("set")) 
    err = lbnode_parse_set(interp, argc-1, argv+1, coord);
    else {
    Tcl_AppendResult(interp, "unknown feature \"", argv[0], "\" of lbnode_extforce", (char *)NULL);
    return  TCL_ERROR;
    }     
  return err;
}

/** Set nemd method to exchange and set nemd parameters */
int tclcommand_nemd_parse_exchange(Tcl_Interp *interp, int argc, char **argv) 
{
  int n_slabs, n_exchange;
  if (argc < 4) {
    Tcl_AppendResult(interp, "wrong # args:  ", (char *)NULL);
    return tclcommand_nemd_print_usage(interp);
  }
  if ( !ARG_IS_I(2, n_slabs) || !ARG_IS_I(3, n_exchange) ) {
    Tcl_AppendResult(interp, "wrong argument type:  ", (char *)NULL);
    return tclcommand_nemd_print_usage(interp);
  }
  /* parameter sanity */
  if ( n_slabs<0 || n_slabs%2!=0 ) {
    Tcl_AppendResult(interp, "nemd <n_slabs> must be non negative and even!",(char *)NULL);
    return (TCL_ERROR);
  }  
  if ( n_slabs > 0 && n_exchange < 0 ) {
    Tcl_AppendResult(interp, "nemd <n_exchange> must be positive!",(char *)NULL);
    return (TCL_ERROR);
  } 

  nemd_method = NEMD_METHOD_EXCHANGE;
  nemd_init(n_slabs, n_exchange, 0.0);
  return (TCL_OK);
}

int tclcommand_localeps(Tcl_Interp* interp, int argc, char** argv)
{
    int mesh = maggs_get_mesh_1D();
    int node_x, node_y, node_z, direction;
    double relative_epsilon;
	
    /* number of arguments has to be 8 */
    if(argc != 9) {
        Tcl_AppendResult(interp, "Wrong number of paramters. Usage: \n", (char *) NULL);
        Tcl_AppendResult(interp, "inter coulomb <bjerrum> memd localeps node <x> <y> <z> dir <X/Y/Z> eps <epsilon>", (char *) NULL);
        return TCL_ERROR;
    }
	
    /* first argument should be "node" */
    if(! ARG_IS_S(1, "node")) return TCL_ERROR;
    
    /* arguments 2-4 should be integers */
    if(! ARG_IS_I(2, node_x)) {
        Tcl_AppendResult(interp, "integer expected", (char *) NULL);
        return TCL_ERROR; }
    if(! ARG_IS_I(3, node_y)) {
        Tcl_AppendResult(interp, "integer expected", (char *) NULL);
        return TCL_ERROR; }
    if(! ARG_IS_I(4, node_z)) {
        Tcl_AppendResult(interp, "integer expected", (char *) NULL);
        return TCL_ERROR; }
    /* check if mesh position is in range */
    if ( (node_x < 0) || (node_y < 0) || (node_z < 0) || (node_x > mesh) || (node_y > mesh) || (node_z > mesh) ) {
        char buffer[TCL_INTEGER_SPACE];
        sprintf(buffer, "%d", mesh);
        Tcl_AppendResult(interp, "epsilon position out of mesh range. Mesh in each dimension is ", buffer, ".", (char *) NULL);
        return TCL_ERROR;
    }
    
    /* parse fifth and sixth argument (e.g. dir X) */
    if(! ARG_IS_S(5, "dir")) return TCL_ERROR;
    if ( (! ARG_IS_S(6, "X")) && (! ARG_IS_S(6, "Y")) && (! ARG_IS_S(6, "Z")) ) {
        Tcl_AppendResult(interp, "Parameter dir should be 'X', 'Y' or 'Z'.", (char *) NULL);
        return TCL_ERROR; }
    if(ARG_IS_S(6, "X")) direction = 0;
    if(ARG_IS_S(6, "Y")) direction = 1;
    if(ARG_IS_S(6, "Z")) direction = 2;
    
    /* parse seventh and eight argument (e.g. eps 0.5) */
    if(! ARG_IS_S(7, "eps")) return TCL_ERROR;
    if ( (! ARG_IS_D(8, relative_epsilon)) || (relative_epsilon < 0.0) ) {
        Tcl_AppendResult(interp, "eps expects a positive double", (char *) NULL);
        return TCL_ERROR; }
    
    double eps_before = maggs_set_permittivity(node_x, node_y, node_z, direction, relative_epsilon);
    
    if (eps_before == 1.0) return TCL_OK;
    else return TCL_OK;
}

static int tclcommand_analyze_fluid_parse_velprof(Tcl_Interp *interp, int argc, char **argv) {
    int i, pdir, vcomp, x1, x2;
    char buffer[TCL_DOUBLE_SPACE];
    double *velprof;

    //fprintf(stderr, "NOTE: analyze fluid velprof is not completely implemented by now.\n      The calling interface might still change without backwards compatibility!\n");

    /*
    if (n_nodes > 1) {
	Tcl_AppendResult(interp, "velocity profile not yet implemented for parallel execution!", (char *)NULL);
	return TCL_ERROR;
    }
    */

    if (argc < 4) {
	Tcl_AppendResult(interp, "usage: analyze fluid velprof <v_comp> <p_dir> <x1> <x2>", (char *)NULL);
	return TCL_ERROR;
    }

    if (!ARG_IS_I(0,vcomp)) return TCL_ERROR;
    if (!ARG_IS_I(1,pdir)) return TCL_ERROR;
    if (!ARG_IS_I(2,x1)) return TCL_ERROR;
    if (!ARG_IS_I(3,x2)) return TCL_ERROR;

    if (pdir != 2) {
	Tcl_AppendResult(interp, "analyze fluid velprof is only implemented for pdir=2 yet!", (char *)NULL);
	return TCL_ERROR;
    }
    if (vcomp != 0) {
	Tcl_AppendResult(interp, "analyze fluid velprof is only implemented for vdir=0 yet", (char *)NULL);
	return TCL_ERROR;
    }

    velprof = malloc(box_l[pdir]/lblattice.agrid*sizeof(double));

    lb_master_calc_velprof(velprof, vcomp, pdir, x1, x2);

    for (i=0; i<box_l[pdir]/lblattice.agrid; i++) {
	Tcl_PrintDouble(interp, i*lblattice.agrid, buffer);
	Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
	Tcl_PrintDouble(interp, velprof[i], buffer);
	Tcl_AppendResult(interp, buffer, "\n", (char *)NULL);
    }	

    free(velprof);

    return TCL_OK;

}

int tclcommand_observable_average(Tcl_Interp* interp, int argc, char** argv, int* change, observable* obs) {
  int reference_observable;
  if (argc < 2) {
    Tcl_AppendResult(interp, "observable new average <reference_id>", (char *)NULL);
    return TCL_ERROR;
  }
  if (!ARG_IS_I(1,reference_observable)) {
    Tcl_AppendResult(interp, "observable new average <reference_id>", (char *)NULL);
    return TCL_ERROR;
  }
  if (reference_observable >= n_observables) {
    Tcl_AppendResult(interp, "The reference observable does not exist.", (char *)NULL);
    return TCL_ERROR;
  }
  observable_average_container* container=(observable_average_container*)malloc(sizeof(observable_average_container));
  container->reference_observable = observables[reference_observable];
  container->n_sweeps = 0;
  obs->n = container->reference_observable->n;
  obs->last_value=(double*)malloc(obs->n*sizeof(double));
  for (int i=0; i<obs->n; i++) 
    obs->last_value[i] = 0;

  obs->container=container;
  obs->update=&observable_update_average;
  obs->calculate=0;

  return TCL_OK;
}

int tclcommand_analyze_set_parse_chain_topology(Tcl_Interp *interp, int argc, char **argv)
{
  /* parses a chain topology (e.g. in 'analyze ( rg | <rg> ) [chain start n chains chain length]' , or
     in 'analyze set chains <chain_start> <n_chains> <chain_length>') */
  int m, i, pc;
  
  if (argc < 3) {
    Tcl_AppendResult(interp, "chain structure info consists of <start> <n> <length>", (char *)NULL);    
    return TCL_ERROR;
  }

  if (! (ARG0_IS_I(chain_start) && ARG1_IS_I(chain_n_chains) && ARG_IS_I(2, chain_length)))
    return TCL_ERROR;

  realloc_topology(chain_n_chains);
  pc = 0;
  for (m = 0; m < n_molecules; m++) {
    topology[m].type = 0;
    realloc_intlist(&topology[m].part, topology[m].part.n = chain_length);
    for (i = 0; i < chain_length; i++)
      topology[m].part.e[i] = pc++;
  }
 
  return TCL_OK;
}

int tclcommand_integrate(ClientData data, Tcl_Interp *interp, int argc, char **argv) 
{
  int  n_steps;
  
  INTEG_TRACE(fprintf(stderr,"%d: integrate:\n",this_node));

  if (argc < 1) {
    Tcl_AppendResult(interp, "wrong # args: \n\"", (char *) NULL);
    return tclcommand_integrate_print_usage(interp);  }
  else if (argc < 2) {                    return tclcommand_integrate_print_status(interp); }

  if (ARG1_IS_S("set")) {
    if      (argc < 3)                    return tclcommand_integrate_print_status(interp);
    if      (ARG_IS_S(2,"nvt"))           return tclcommand_integrate_set_nvt(interp, argc, argv);
#ifdef NPT
    else if (ARG_IS_S(2,"npt_isotropic")) return tclcommand_integrate_set_npt_isotropic(interp, argc, argv);
#endif
    else {
      Tcl_AppendResult(interp, "unknown integrator method:\n", (char *)NULL);
      return tclcommand_integrate_print_usage(interp);
    }
  } else if ( !ARG_IS_I(1,n_steps) ) return tclcommand_integrate_print_usage(interp);

  /* go on with integrate <n_steps> */
  if(n_steps < 0) {
    Tcl_AppendResult(interp, "illegal number of steps (must be >0) \n", (char *) NULL);
    return tclcommand_integrate_print_usage(interp);;
  }
  /* perform integration */
  if (mpi_integrate(n_steps))
    return mpi_gather_runtime_errors(interp, TCL_OK);
  return TCL_OK;
}

int tclcommand_inter_parse_comfixed(Tcl_Interp * interp,
				    int part_type_a, int part_type_b,
				    int argc, char ** argv)
{
  int flagc;
  	
  if (argc != 2) {
    Tcl_AppendResult(interp, "comfixed needs 1 parameters: "
		     "<comfixed_flag> ", (char *) NULL);
    return 0;
  }
	 
  if (part_type_a != part_type_b) {
    Tcl_AppendResult(interp, "comfixed must be among same type interactions", (char *) NULL);
    return 0;
  }

  /* copy comfixed parameters */
  if ((! ARG_IS_I(1, flagc)) ) {
    Tcl_AppendResult(interp, "comfixed needs 1 INTEGER parameter: "
		     "<comfixed_flag>", (char *) NULL);
    return 0;
  }

  switch (comfixed_set_params(part_type_a, part_type_b, flagc)) {
  case 1:
    Tcl_AppendResult(interp, "particle types must be non-negative", (char *) NULL);
    return 0;
  case 2:
    Tcl_AppendResult(interp, "works only with a single CPU", (char *) NULL);
    return 0;
  }

  return 2;
}

int tclcommand_thermostat_parse_inter_dpd(Tcl_Interp *interp, int argc, char ** argv)
{
  double temp;

  if (argc < 2) {
    Tcl_AppendResult(interp, "thermostat needs 1 parameter: "
		     "<temperature>",
		     (char *) NULL);
    return TCL_ERROR;
  }

  if (argc>2 && ARG_IS_S(2, "ignore_fixed_particles")) {
    if (argc == 3)
      dpd_ignore_fixed_particles=1;
    else if (argc!= 4 || (!ARG_IS_I(3, dpd_ignore_fixed_particles))) 
      return TCL_ERROR;
    mpi_bcast_parameter(FIELD_DPD_IGNORE_FIXED_PARTICLES);
    return TCL_OK;
  }

  /* copy lattice-boltzmann parameters */
  if (! ARG_IS_D(2, temp)) { return TCL_ERROR; }

  if ( temp < 0.0 ) {
    Tcl_AppendResult(interp, "temperature must be non-negative", (char *) NULL);
    return TCL_ERROR;
  }
  temperature = temp;
  thermo_switch = ( thermo_switch | THERMO_INTER_DPD );
  mpi_bcast_parameter(FIELD_THERMO_SWITCH);
  mpi_bcast_parameter(FIELD_TEMPERATURE);
  return (TCL_OK);
}

int tclcommand_inter_parse_interrf(Tcl_Interp * interp,
				   int part_type_a, int part_type_b,
				   int argc, char ** argv)
{
  /* parameters needed for RF */
  int rf_on;
  int change;

  /* get reaction_field interaction type */
  if (argc < 2) {
    Tcl_AppendResult(interp, "inter_rf needs 1 parameter: "
		     "<rf_on>",
		     (char *) NULL);
    return 0;
  }

  /* copy reaction_field parameters */
  if (! ARG_IS_I(1, rf_on)) {
    Tcl_AppendResult(interp, "<rf_on> must be int",
		     (char *) NULL);
    return 0;
  }
  change = 2;
	
  if (! ((rf_on==0) || (rf_on==1)) ) {
    Tcl_AppendResult(interp, "rf_on must be 0 or 1", (char *) NULL);
    return 0;
  }
  if (interrf_set_params(part_type_a, part_type_b,rf_on) == ES_ERROR) {
    Tcl_AppendResult(interp, "particle types must be non-negative", (char *) NULL);
    return 0;
  }
  return change;
}

int tclcommand_thermostat_parse_cpu(Tcl_Interp *interp, int argc, char **argv) 
{
  int temp;

  #ifndef __linux__
  Tcl_AppendResult(interp, "This feature is currently only supported on Linux platforms.", (char *)NULL);
  return (TCL_ERROR);
  #endif

  /* check number of arguments */
  if (argc < 3) {
    Tcl_AppendResult(interp, "wrong # args:  should be \n\"",
         argv[0]," ",argv[1]," <temp>\"", (char *)NULL);
    return (TCL_ERROR);
  }

  /* check argument types */
  if ( !ARG_IS_I(2, temp) ) {
    Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs one INT", (char *)NULL);
    return (TCL_ERROR);
  }

  /* broadcast parameters */
  temperature = temp;
  thermo_switch = ( thermo_switch | THERMO_CPU );
  mpi_bcast_parameter(FIELD_THERMO_SWITCH);
  mpi_bcast_parameter(FIELD_TEMPERATURE);
  return (TCL_OK);
}

int tclcommand_minimize_energy(ClientData data, Tcl_Interp *interp, int argc, char **argv) 
{
  int  max_steps;
  double f_max, gamma, max_displacement;
  
  if (argc != 5) {
    Tcl_AppendResult(interp, "wrong # args: \n\"", (char *) NULL);
    return usage(interp);
  }
  else {
    if(!ARG_IS_D(1,f_max)) {
      return usage(interp);      
    }
    if(!ARG_IS_I(2,max_steps)) {
      return usage(interp);      
    }
    if(!ARG_IS_D(3,gamma)) {
      return usage(interp);      
    }
    if(!ARG_IS_D(4,max_displacement)) {
      return usage(interp);      
    }
  }

  minimize_energy_init(f_max, gamma, max_steps, max_displacement);
  mpi_minimize_energy();

  return TCL_OK;
}

int tclcommand_inter_parse_ljgen(Tcl_Interp * interp,
                                 int part_type_a, int part_type_b,
                                 int argc, char ** argv)
{
    /* parameters needed for LJGEN */
    double eps, sig, cut, shift, offset, cap_radius, b1, b2;
    int change, a1, a2;

    /* get lennard-jones interaction type */
    if (argc < 10) {
        Tcl_AppendResult(interp, "lj-gen needs 9 parameters: "
                         "<lj_eps> <lj_sig> <lj_cut> <lj_shift> <lj_offset> <a1> <a2> <b1> <b2>",
                         (char *) NULL);
        return 0;
    }

    /* copy lennard-jones parameters */
    if ((! ARG_IS_D(1, eps))    ||
            (! ARG_IS_D(2, sig))    ||
            (! ARG_IS_D(3, cut))    ||
            (! ARG_IS_D(4, shift))  ||
            (! ARG_IS_D(5, offset)) ||
            (! ARG_IS_I(6, a1))     ||
            (! ARG_IS_I(7, a2))     ||
            (! ARG_IS_D(8, b1))     ||
            (! ARG_IS_D(9, b2))) {
        Tcl_AppendResult(interp, "lj-gen needs 7 DOUBLE and 2 INT parameers: "
                         "<lj_eps> <lj_sig> <lj_cut> <lj_shift> <lj_offset> <a1> <a2> <b1> <b2>",
                         (char *) NULL);
        return ES_ERROR;
    }
    change = 10;

    cap_radius = -1.0;
    /* check wether there is an additional double, cap radius, and parse in */
    if (argc >= 11 && ARG_IS_D(10, cap_radius))
        change++;
    else
        Tcl_ResetResult(interp);
    if (ljgen_set_params(part_type_a, part_type_b,
                         eps, sig, cut, shift, offset, a1, a2, b1, b2,
                         cap_radius) == ES_ERROR) {
        Tcl_AppendResult(interp, "particle types must be non-negative", (char *) NULL);
        return 0;
    }
    return change;
}

int tclcommand_metadynamics_parse_relative_z(Tcl_Interp *interp, int argc, char **argv)
{
  int    pid1, pid2, dbins, numrelaxationsteps;
  double dmin, dmax, bheight, bwidth, fbound;

  /* check number of arguments */
  if (argc < 11) {
    Tcl_AppendResult(interp, "wrong # args:  should be \n\"",
                     argv[0]," ",argv[1]," <pid1> <pid2> <z_min> <z_max> <b_height> <b_width> <f_bound> <z_bins> <num_relaxation_steps>\"", (char *)NULL);
    return (TCL_ERROR);
  }

  /* check argument types */
  if ( !ARG_IS_I(2, pid1) || !ARG_IS_I(3, pid2) || !ARG_IS_D(4, dmin) || !ARG_IS_D(5, dmax) ||
       !ARG_IS_D(6, bheight) || !ARG_IS_D(7, bwidth) || !ARG_IS_D(8, fbound) || !ARG_IS_I(9, dbins) || !ARG_IS_I(10, numrelaxationsteps) ) {
    Tcl_AppendResult(interp, argv[0]," ",argv[1]," needs two INTS, five DOUBLES, and two INTS in this order", (char *)NULL);
    return (TCL_ERROR);
  }

  if (pid1 < 0 || pid1 > max_seen_particle || pid2 < 0 || pid2 > max_seen_particle) {
    Tcl_AppendResult(interp, "pid1 and/or pid2 out of range", (char *)NULL);
    return (TCL_ERROR);
  }

  if (dmax < dmin || bheight < 0 || bwidth < 0 || fbound < 0 || dbins < 0 || numrelaxationsteps <0) {
    Tcl_AppendResult(interp, "check parameters: inconcistency somewhere", (char *)NULL);
    return (TCL_ERROR);
  }

  free(meta_acc_force);
  free(meta_acc_fprofile);

  /* broadcast parameters */
  meta_pid1         = pid1;
  meta_pid2         = pid2;
  meta_bias_height  = bheight;
  meta_bias_width   = bwidth;
  meta_xi_min       = dmin;
  meta_xi_max       = dmax;
  meta_f_bound      = fbound;
  meta_xi_num_bins  = dbins;

  meta_switch = META_REL_Z;
  meta_num_relaxation_steps = numrelaxationsteps;

  return (TCL_OK);
}

int tclcommand_inter_coulomb_parse_ewald(Tcl_Interp * interp, int argc, char ** argv)
{
  double r_cut, alpha;
  int i, kmax;

  coulomb.method = COULOMB_EWALD;
    
#ifdef PARTIAL_PERIODIC
  if(PERIODIC(0) == 0 ||
     PERIODIC(1) == 0 ||
     PERIODIC(2) == 0)
    {
      Tcl_AppendResult(interp, "Need periodicity (1,1,1) with Coulomb EWALD",
		       (char *) NULL);
      return TCL_ERROR;  
    }
#endif

  if (argc < 2) {
    Tcl_AppendResult(interp, "expected: inter coulomb <bjerrum> ewald <r_cut> <alpha> <kmax>",
		     (char *) NULL);
    return TCL_ERROR;  
  }

  if(! ARG0_IS_D(r_cut))
    return TCL_ERROR;  

  if(argc != 3) {
    Tcl_AppendResult(interp, "wrong # arguments: inter coulomb <bjerrum> ewald <r_cut> <alpha> <kmax>",
		     (char *) NULL);
    return TCL_ERROR;  
  }

  if(! ARG_IS_D(1, alpha))
    return TCL_ERROR;

  if(! ARG_IS_I(2, kmax))
    return TCL_ERROR;

  if ((i = ewald_set_params(r_cut, alpha, kmax)) < 0) {
    switch (i) {
    case -1:
      Tcl_AppendResult(interp, "r_cut must be positive", (char *) NULL);
      break;
    case -4:
      Tcl_AppendResult(interp, "alpha must be positive", (char *) NULL);
      break;
    case -5:
      Tcl_AppendResult(interp, "kmax must be greater than zero", (char *) NULL);
    default:;
      Tcl_AppendResult(interp, "unspecified error", (char *) NULL);
    }

    return TCL_ERROR;

  }

  return TCL_OK;
}

/** parse TCL command.
    number of parameters is checked and maggs_set_parameters function is called.
    @return zero if successful
    @param interp  TCL interpreter handle
    @param argc    number of arguments given
    @param argv    array of arguments given
*/
int tclcommand_inter_coulomb_parse_maggs(Tcl_Interp * interp, int argc, char ** argv)
{
    int mesh;
    double f_mass;
    double epsilon = 1.0;
    int finite_epsilon_flag = 1;
	
    /* if the command is localeps, call function */
    if ( (argc > 0) && (ARG_IS_S(0, "localeps")) )
        return tclcommand_localeps(interp, argc, argv);
    
    if(argc < 2) {
        Tcl_AppendResult(interp, "Not enough parameters: inter coulomb <bjerrum> memd <f_mass> <mesh>", (char *) NULL);
        return TCL_ERROR;
    }
	
    if(! ARG_IS_D(0, f_mass))
        return TCL_ERROR;
	
    if(! ARG_IS_I(1, mesh)) {
        Tcl_AppendResult(interp, "integer expected", (char *) NULL);
        return TCL_ERROR;
    }
	
    if(argc > 4) {
        Tcl_AppendResult(interp, "Too many parameters: inter coulomb memd <f_mass> <mesh> [epsilon <eps>]", (char *) NULL);
        return TCL_ERROR;
    }
    if(argc == 3) {
        Tcl_AppendResult(interp, "Usage: inter coulomb memd <f_mass> <mesh> [epsilon <eps>]", (char *) NULL);
        return TCL_ERROR;
    }
    if(argc == 4) {
        if (ARG_IS_S(2, "epsilon")) {
            if(! (ARG_IS_D(3, epsilon) && epsilon > 0.0)) {
                Tcl_AppendResult(interp, "epsilon expects a positive double",
                             (char *) NULL);
                return TCL_ERROR;
            }
        }
    } else finite_epsilon_flag=1;

  coulomb.method = COULOMB_MAGGS;
	
  int res = maggs_set_parameters(coulomb.bjerrum, f_mass, mesh,
                                 finite_epsilon_flag, epsilon);
  switch (res) {
  case -1:
    Tcl_AppendResult(interp, "mass of the field is negative", (char *)NULL);
    return TCL_ERROR;
  case -2:
    Tcl_AppendResult(interp, "mesh must be positive", (char *) NULL);
    return TCL_ERROR;
  case ES_OK:
    return TCL_OK;
  }
  Tcl_AppendResult(interp, "unknown error", (char *) NULL);
  return TCL_ERROR;
}

int tclcommand_cellsystem(ClientData data, Tcl_Interp *interp,
	       int argc, char **argv)
{
  int err = 0;

  if (argc <= 1) {
    Tcl_AppendResult(interp, "usage: cellsystem <system> <params>", (char *)NULL);
    return TCL_ERROR;
  }

  if (ARG1_IS_S("domain_decomposition")) {
    if (argc > 2) {
      if (ARG_IS_S(2,"-verlet_list"))
	dd.use_vList = 1;
      else if(ARG_IS_S(2,"-no_verlet_list")) 
	dd.use_vList = 0;
      else{
	Tcl_AppendResult(interp, "wrong flag to",argv[0],
			 " : should be \" -verlet_list or -no_verlet_list \"",
			 (char *) NULL);
	return (TCL_ERROR);
      }
    }
    /** by default use verlet list */
    else dd.use_vList = 1;
    mpi_bcast_cell_structure(CELL_STRUCTURE_DOMDEC);
  }
  else if (ARG1_IS_S("nsquare"))
    mpi_bcast_cell_structure(CELL_STRUCTURE_NSQUARE);
  else if (ARG1_IS_S("layered")) {
    if (argc > 2) {
      if (!ARG_IS_I(2, n_layers))
	return TCL_ERROR;
      if (n_layers <= 0) {
	Tcl_AppendResult(interp, "layer height should be positive", (char *)NULL);
	return TCL_ERROR;
      }
      determine_n_layers = 0;
    }

    /* check node grid. All we can do is 1x1xn. */
    if (node_grid[0] != 1 || node_grid[1] != 1) {
      node_grid[0] = node_grid[1] = 1;
      node_grid[2] = n_nodes;
      
      err = mpi_bcast_parameter(FIELD_NODEGRID);
    }
    else
      err = 0;

    if (!err)
      mpi_bcast_cell_structure(CELL_STRUCTURE_LAYERED);
  }
  else {
    Tcl_AppendResult(interp, "unkown cell structure type \"", argv[1],"\"", (char *)NULL);
    return TCL_ERROR;
  }
  return gather_runtime_errors(interp, TCL_OK);
}

int tclcommand_inter_parse_inter_dpd(Tcl_Interp * interp,
		       int part_type_a, int part_type_b,
		       int argc, char ** argv)
{
  /* parameters needed for LJ */
  extern double temperature;
  double gamma,r_c,tgamma,tr_c;
  int wf,twf;
  int change;

  /* get inter_dpd interaction type */
  if (argc < 7) {
    Tcl_AppendResult(interp, "inter_dpd needs 6 parameters: "
		     "<gamma> <r_cut> <wf> <tgamma> <tr_cut> <twf>",
		     (char *) NULL);
    return 0;
  }
  if (temperature == -1) {
    Tcl_AppendResult(interp, "Please set temperature first:  temperature inter_dpd temp",(char *) NULL);
    return 0;
  }

  /* copy lennard-jones parameters */
  if ((! ARG_IS_D(1, gamma))  ||
      (! ARG_IS_D(2, r_c))    ||
      (! ARG_IS_I(3, wf))     ||
      (! ARG_IS_D(4, tgamma)) ||
      (! ARG_IS_D(5, tr_c))    ||
      (! ARG_IS_I(6, twf))    ) {
    Tcl_AppendResult(interp, "inter_dpd needs 6 parameters: "
		     "<gamma> <r_cut> <wf> <tgamma> <tr_cut> <twf> ",
		     (char *) NULL);
    return 0;
  }
  change = 7;
	
  if (inter_dpd_set_params(part_type_a, part_type_b,
			       gamma,r_c,wf,tgamma,tr_c,twf) == ES_ERROR) {
    Tcl_AppendResult(interp, "particle types must be non-negative", (char *) NULL);
    return 0;
  }
  inter_dpd_init();

  return change;

}

int tclcommand_analyze_parse_formfactor(Tcl_Interp *interp, int average, int argc, char **argv)
{
  /* 'analyze { formfactor | <formfactor> } <qmin> <qmax> <qbins> [<chain_start> <n_chains> <chain_length>]' */
  /***********************************************************************************************************/
  char buffer[2*TCL_DOUBLE_SPACE+5];
  int i;
  double qmin,qmax, q,qfak, *ff; int qbins;
  if (argc < 3) {
    Tcl_AppendResult(interp, "Wrong # of args! Usage: analyze formfactor <qmin> <qmax> <qbins> [<chain_start> <n_chains> <chain_length>]",
		     (char *)NULL);
    return (TCL_ERROR);
  } else {
    if (!ARG0_IS_D(qmin))
      return (TCL_ERROR);
    if (!ARG1_IS_D(qmax))
      return (TCL_ERROR);
    if (!ARG_IS_I(2, qbins))
      return (TCL_ERROR);
    argc-=3; argv+=3;
  }
  if (tclcommand_analyze_set_parse_chain_topology_check(interp, argc, argv) == TCL_ERROR) return TCL_ERROR;

  if ((chain_n_chains == 0) || (chain_length == 0)) {
    Tcl_AppendResult(interp, "The chain topology has not been set",(char *)NULL); return TCL_ERROR;
  }
  
  if (qbins <=0) {
    Tcl_AppendResult(interp, "Nothing to be done - choose <qbins> greater zero to get S(q)!",(char *)NULL); return TCL_ERROR;
  }

  if (qmin <= 0.) {
    Tcl_AppendResult(interp, "formfactor S(q) requires qmin > 0", (char *)NULL);
    return TCL_ERROR;
  }
  if (qmax <= qmin) {
    Tcl_AppendResult(interp, "formfactor S(q) requires qmin < qmax", (char *)NULL);
    return TCL_ERROR;
  }

  if (!average) analyze_formfactor(qmin, qmax, qbins, &ff);
  else if (n_configs == 0) {
    Tcl_AppendResult(interp, "no configurations found! ", (char *)NULL);
    Tcl_AppendResult(interp, "Use 'analyze append' to save some, or 'analyze formfactor ...' to only look at current state!",
		     (char *)NULL);
    return TCL_ERROR; }
  else analyze_formfactor_av(qmin, qmax, qbins, &ff);
  
  q = qmin; qfak = pow((qmax/qmin),(1.0/qbins));
  for(i=0; i<=qbins; i++) { sprintf(buffer,"{%f %f} ",q,ff[i]); q*=qfak; Tcl_AppendResult(interp, buffer, (char *)NULL); }
  free(ff);
  return (TCL_OK);
}