Re: [ESPResSo] Position-dependent Langevoin thermostat

[Lorenzo Isella]

>> I would like to create somewhere two text files to be read (only once)
>> by Espresso into two arrays eta1[] and eta[2] use them as explained above.
>> Maybe I should add that I aim at having e.g. two (or more) sets of 50
>> monomers each (each set would be an aggregate; let us call them  A and
>> B) and having eta1[] and eta2[]  as arrays of length 50. Then I would
>> use the same eta1(2) correction for the monomers in A and B since the
>> two aggregates are identical and their corresponding monomers have the
>> same number of first neighbors.

>>
>
> If you choose eta1 and eta2 for each particle, things are relatively easy.
> The best way to implement this is to add a new particle property "etas";
> take a look at e.g. how the charge is implemented, for example. First, you
> add an array double eta[2] to ParticleProperties in particle_data.h. Then,
> you need to write two new routines in particle_data.c, part_print_etas and
> part_parse_etas, to set and read the values from Tcl. You do this similar to
> the same ones for the charge, just that you read/write two doubles. Finally,
> you need to initialize your values in init_particle in particle_data.c.
> After that, you can set the values from Tcl. So, you open the text files in
> Tcl and read them line by line into the particles.
>
>
> Many regards,
> Axel
>

Dear Axel,
Thanks a lot for your support! I am now implementing the modifications
you suggested.
I think I am comfortable with most of them (or at least I need to test
them before posting again), but there is something that puzzles me.
When I look at the routine to parse q (which I paste in the following)
that is the "model" routine I should modify

int part_parse_q(Tcl_Interp *interp, int argc, char **argv,
                 int part_num, int * change)
{
    double q;

    *change = 1;

    if (argc < 1) {
      Tcl_AppendResult(interp, "q requires 1 argument", (char *) NULL);
      return TCL_ERROR;
    }

    /* set charge */
    if (! ARG0_IS_D(q))
      return TCL_ERROR;

    if (set_particle_q(part_num, q) == TCL_ERROR) {
      Tcl_AppendResult(interp, "set particle position first", (char *)NULL);

      return TCL_ERROR;
    }

    return TCL_OK;
}



I see that there is a call to a set_particle_q function

int set_particle_q(int part, double q)
{
  int pnode;
  if (!particle_node)
    build_particle_node();

  if (part < 0 || part > max_seen_particle)
    return TCL_ERROR;
  pnode = particle_node[part];

  if (pnode == -1)
    return TCL_ERROR;
  mpi_send_q(pnode, part, q);
  return TCL_OK;
}

which in turn deals with the mpi_send_q function.
I tracked that function as ended up in section of the code dealing
with the inner workings of mpi and parallelization (which I would
definitely like to leave untouched at this stage).
Do I really need to do browse though all this or can a simpler
modification of part_parse_q do what I need?
Many thanks

Lorenzo Isella