Re: [ESPResSo] Position-dependent Langevoin thermostat

[Lorenzo Isella]

Hello,
Sorry posting again on this, but I am completely puzzled at what I am seeing.
I anticipate my problem: if I initialize eta[0]=eta[1]=1.0 in Langevin
thermostat and I do not read anything into eta[] in a simulation, then
nothing can change with respect to using  the standard Langevin
thermostat.
Instead I must be doing something wrong. When I re-ran an old test
simulation of mine,  I simply did not see any "evolution" of the
system any longer (!), but everything stays as in the initial state. I
am surprised because I did not touch the Verlet integrator at all.
BTW, this does not happen of course when I use espresso 2.1 (but I do
not believe that the specific espresso version is the point here) as
downloaded from the espresso web site without any modifications. I am
not posting the simulation + input data unless this turns out to be
necessary.
This is the content of my configuration file myconfig.h

/** \file myconfig-sample.h

    This is a sample for the file myconfig.h.

    Uncomment any of the following lines to myconfig.h to activate
    the corresponding feature of Espresso. It is recommended to turn
    only those features on that you actually need to optimize the
    performance of Espresso for your problem. For details on these
    features see the user's guide, Sec. 2.6 "Configuration options".

    To access the information on the compilation status of the code
    you are working with in your Espresso Tcl-script, use the
    corresponding \ref tcl_features "Tcl-commands".

    If you add a new feature to Espresso, you also have to add the
    corresponding lines in the function \ref compilation_callback and
    to add documentation in <tt>doc/text/features.doc</tt>.

    <b>Responsible:</b>
    <a href="mailto:address@hidden";>Axel</a>

*/

/**********************************************************************/
/*                       general core features                        */
/**********************************************************************/

 #define PARTIAL_PERIODIC
#define ELECTROSTATICS
/* #define ROTATION  */
 #define EXTERNAL_FORCES
 #define CONSTRAINTS
 #define MASS
 #define EXCLUSIONS
 #define COMFORCE
 #define COMFIXED
 #define MOLFORCES
 #define BOND_CONSTRAINT
 #define MODES

/**********************************************************************/
/*                        integrator features                         */
/**********************************************************************/

 #define NEMD
 #define NPT
 #define DPD

/**********************************************************************/
/*                           interactions                             */
/**********************************************************************/

 #define TABULATED
 #define LENNARD_JONES
 #define SMOOTH_STEP
 #define BMHTF_NACL
#define LJ_WARN_WHEN_CLOSE
 #define MORSE
 #define LJCOS
 #define LJCOS2
 #define BUCKINGHAM
 #define SOFT_SPHERE

/* Note: Activate ONLY ONE bonded angle potential out of the following! */
/* #define BOND_ANGLE_HARMONIC */
 #define BOND_ANGLE_COSINE
/* #define BOND_ANGLE_COSSQUARE */

/**********************************************************************/
/*                            debugging                               */
/**********************************************************************/

/* #define ADDITIONAL_CHECKS

 #define COMM_DEBUG
 #define EVENT_DEBUG
#define INTEG_DEBUG
 #define CELL_DEBUG
 #define GHOST_DEBUG
 #define LATTICE_DEBUG
 #define HALO_DEBUG
 #define GRID_DEBUG
 #define VERLET_DEBUG
 #define PARTICLE_DEBUG
 #define P3M_DEBUG
 #define EWALD_DEBUG
 #define FFT_DEBUG
 #define RANDOM_DEBUG
 #define FORCE_DEBUG
 #define THERMO_DEBUG
 #define LJ_DEBUG
 #define MORSE_DEBUG */

/* #define ESR_DEBUG */
/* #define ESK_DEBUG */
/* #define FENE_DEBUG */
/* #define GHOST_FORCE_DEBUG */
/* #define ONEPART_DEBUG 13 */
/* #define STAT_DEBUG */
/* #define POLY_DEBUG */
/* #define MOLFORCES_DEBUG */
/* #define MEM_DEBUG */
/* #define MAGGS_DEBUG */

/* #define ASYNC_BARRIER */

/* #define MPI_CORE */
/* #define FORCE_CORE */

>From the directory where I am building espresso 2.1 with my
modifications I issue this command (containing the path to the
modified espresso source code)

$ /media/disk/espresso_modified/espresso-desperation/configure
CPPFLAGS=-I/usr/include/tcl8.4 --without-mpi

followed by

$ make -j 4

I now post the output of a few diff's to show what modifications I
actually implemented



diff  /media/disk/espresso_modified/espresso-desperation/thermostat.c
/media/disk/espresso_modified/espresso-2.1.0/thermostat.c

< /* Add the eta array */
<
< double eta[2] ;
<
<
270,271c265,266
<   langevin_pref1 = -eta[0]*langevin_gamma/time_step;
<   langevin_pref2 = sqrt(eta[1]*24.0*temperature*langevin_gamma/time_step);
---
>   langevin_pref1 = -langevin_gamma/time_step;
>   langevin_pref2 = sqrt(24.0*temperature*langevin_gamma/time_step);

I did not touch thermostat.h, while I did some work in particle data.c
(one routine to print and one parse eta[]  and one to set it without
resorting to  MPI, as you suggested)

diff  /media/disk/espresso_modified/espresso-desperation/particle_data.c
/media/disk/espresso_modified/espresso-2.1.0/particle_data.c
92,97d91
<   part->p.eta[0]   = 1.0;
<   part->p.eta[1]   = 1.0;
<
<
<
<
482,499d475
<
<
< /* I add a routine to deal with the etas */
<
< void part_print_etas(Particle *part, char *buffer, Tcl_Interp *interp)
< {
<
<   Tcl_PrintDouble(interp, part->p.eta[0], buffer);
<   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
<   Tcl_PrintDouble(interp, part->p.eta[1], buffer);
<   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
< }
<
<
<
<
<
<
1186,1230d1161
<
<
<
<
< int part_parse_etas(Tcl_Interp *interp, int argc, char **argv,
<                int part_num, int * change)
< {
<   double eta[2];
<
<   *change = 2;
<
<   if (argc < 2) {
<     Tcl_AppendResult(interp, "eta requires 2 arguments", (char *) NULL);
<
<     return TCL_ERROR;
<   }
<
<   /* set eta */
<   if (! ARG_IS_D(0, eta[0]))
<     return TCL_ERROR;
<
<   if (! ARG_IS_D(1, eta[1]))
<     return TCL_ERROR;
<
<
<
<   if (set_particle_etas(part_num, eta) == TCL_ERROR) {
<     Tcl_AppendResult(interp, "set particle position first", (char *)NULL);
<
<     return TCL_ERROR;
<   }
<
<   return TCL_OK;
< }
<
<
<
<
<
<
<
<
<
<
<
1983,2009d1913
<
< int set_particle_etas(int part, double eta[2])
< {
<   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_v(pnode, part, eta); */
<
<   Particle *p = local_particles[part];
<
<   p->p.eta[0]=eta[0] ;
<   p->p.eta[1]=eta[1] ;
<
<
<
<   return TCL_OK;
< }
<
<
<

Finally, I added an array and function declaration in particle_data.h


diff  /media/disk/espresso_modified/espresso-desperation/particle_data.h
/media/disk/espresso_modified/espresso-2.1.0/particle_data.h
88,92d87
<
<   /** add an array which contains eta1 and eta */
<
<   double eta[2] ;
<
390,397d384
<
<
< // I added a function declaration to deal with eta
<
< int set_particle_etas(int part, double eta[2]);
<
<


As I wrote above, the problem is that, given the same configuration
file and the same procedure to build espresso, if I use the build from
the modified source code and from 2.1.0 to perform the the same
simulation (which never touches the values in eta) I get two different
results.
Since in the former case I see no motion at all (and I am simulating
the motion of clusters), I know for a fact that my modifications must
have introduced a bug somewhere in espresso, but I do not understand
where and why.
Any suggestion here is really helpful.
Cheers

Lorenzo




2008/11/5 Axel Arnold <address@hidden>:
>>
>> 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)
>>
>> I see that there is a call to a set_particle_q function
>>
>> int set_particle_q(int part, double q)
>>
>> 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?
>
> Ooops, sorry, yes, I somehow forgot that bit, ahem. Now, things depend on
> whether you want to run on more than one CPU or not. If not, you can sort
> of bypass MPI, but things will fail terribly if you ever start your program
> with more than one processor; however, if you compile with --without-mpi,
> that cannot happen.
>
> In this case, you can simply replace the mpi_send_q with
>
> Particle *p = local_particles[part]; p->p.eta1 = eta1;
>
> or so in set_particle_q. In any case, you need to also copy set_particle_q.
>
> If you want to run in parallel, that is in fact not much more difficult.
> You just need to continue copy&paste for set_particle_q, mpi_send_q and
> mpi_send_q_slave, and also add the declaration of mpi_send_q to to
> communication.h.
>
> However, you also need to register the mpi_send_eta_slave
> function (which is invoked on remote processors to retrieve the updated
> eta of a particle). That is done in communication.c: First, you need a new
> serial number for mpi_send_eta. This are all the defines at the beginning.
> Look for #define REQ_MAXIMUM, and assign this value for your code, using
> REQ_ETA for example; increase REQ_MAXIMUM by one, of course.
> Then, add a declaration for mpi_send_eta_slave to the list right below
> the serial numbers, and to the slave_callbacks array. There it needs to be
> at the place indicated by its serial number, that is, at the end.
> Finally, add a fancy name to the names[] array right below, and finished.
> Now, your value can be set for any particle on as many nodes as you wish
> :-).
>
> Cheers,
> Axel
>
> --
> Dr. Axel Arnold
> Fraunhofer SCAI
> Schloss Birlinghoven, 53754 Sankt Augustin, Germany
> Tel: +49 2241 14 2575
>



-- 
Free will does not consist in inverting the river flow, but in being
the fish that leaps upstream.