From f86fbdec0b0c59e52a769be89d4c24715bcbd581 Mon Sep 17 00:00:00 2001 From: "J. de Graaf" Date: Mon, 3 Dec 2012 10:53:03 +0100 Subject: [PATCH 1/5] testing --- scripts/auxiliary.tcl | 14 ++-- src/Makefile.am | 2 + src/initialize.c | 1 + src/particle_data.c | 4 ++ src/particle_data.h | 4 ++ src/reaction.c | 66 +++++++++++++++--- src/reaction.h | 1 + src/tcl/initialize_interpreter.c | 8 +++ src/tcl/reaction_tcl.c | 116 +++++++++++++++++-------------- 9 files changed, 148 insertions(+), 68 deletions(-) diff --git a/scripts/auxiliary.tcl b/scripts/auxiliary.tcl index 02f245799ce..e4a5ade9bdc 100644 --- a/scripts/auxiliary.tcl +++ b/scripts/auxiliary.tcl @@ -474,13 +474,13 @@ proc system_com { } { set npart [setmd n_part] #calculate center of mass - set com {0 0 0} + set com {0.0 0.0 0.0} set part_cnt 0 if {[has_feature MASS]} { set net_mass 0 - for {set i 0} {$part_cnt < $npart} {incr i} { + for {set i 0} {$i < $npart} {incr i} { if {[part $i] != "na"} { set pos [part $i print p] set mass [part $i print mass] @@ -489,10 +489,10 @@ proc system_com { } { incr part_cnt } } - + return [vecscale [expr 1.0/$net_mass] $com] } else { - for {set i 0} {$part_cnt < $npart} {incr i} { + for {set i 0} {$i < $npart} {incr i} { if {[part $i] != "na"} { set pos [part $i print p] set com [vecadd $com $pos] @@ -522,7 +522,7 @@ proc system_com_vel {} { if {[has_feature MASS]} { set net_mass 0 - for {set i 0} {$part_cnt < $npart } {incr i} { + for {set i 0} {$i < $npart} {incr i} { if { [part $i] != "na" } { set part_vel [part $i print v] set part_mass [part $i print mass] @@ -532,9 +532,9 @@ proc system_com_vel {} { } } - return [vecscale [expr 1.0/$net_mass] $com_vel] + return "[vecscale [expr 1.0/$net_mass] $com_vel]" } else { - for {set i 0} {$part_cnt < $npart } {incr i} { + for {set i 0} {$i < $npart} {incr i} { if { [part $i] != "na" } { set part_vel [part $i print v] set com_vel [vecadd $com_vel $part_vel] diff --git a/src/Makefile.am b/src/Makefile.am index 08a7a5e6429..ef0dfdb244f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -45,6 +45,7 @@ libEspresso_a_SOURCES = \ fft-dipolar.c fft-dipolar.h \ forcecap.c forcecap.h \ forces.c forces.h \ + galilei.c galilei.h \ ghosts.c ghosts.h \ global.c global.h \ grid.c grid.h \ @@ -161,6 +162,7 @@ libEspresso_a_SOURCES += \ tcl/constraint_tcl.c tcl/constraint_tcl.h \ tcl/domain_decomposition_tcl.c tcl/domain_decomposition_tcl.h \ tcl/energy_tcl.c \ + tcl/galilei_tcl.c tcl/galilei_tcl.h \ tcl/global_tcl.c tcl/global_tcl.h \ tcl/grid_tcl.c tcl/grid_tcl.h \ tcl/iccp3m_tcl.c tcl/iccp3m_tcl.h \ diff --git a/src/initialize.c b/src/initialize.c index 31d30c74898..c0410dbff49 100644 --- a/src/initialize.c +++ b/src/initialize.c @@ -135,6 +135,7 @@ void on_program_start() #ifdef REACTIONS reaction.back_rate=-1.0; + reaction.sing_mult=0; #endif /* diff --git a/src/particle_data.c b/src/particle_data.c index 0d2c48f827f..69291134171 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -120,6 +120,10 @@ void init_particle(Particle *part) part->p.mu_E[2] = 0.0; #endif +#ifdef REACTIONS + part->p.reacted = 0; +#endif + /* ParticlePosition */ part->r.p[0] = 0.0; part->r.p[1] = 0.0; diff --git a/src/particle_data.h b/src/particle_data.h index 0ad21c21191..ab98a6f6b80 100644 --- a/src/particle_data.h +++ b/src/particle_data.h @@ -124,6 +124,10 @@ typedef struct { double T; double gamma; #endif + +#ifdef REACTIONS + int reacted; +#endif } ParticleProperties; /** Positional information on a particle. Information that is diff --git a/src/reaction.c b/src/reaction.c index cb389faec2a..2c800bc6b16 100644 --- a/src/reaction.c +++ b/src/reaction.c @@ -36,6 +36,7 @@ void setup_reaction() { MPI_Bcast(&reaction.range, 1, MPI_DOUBLE, 0, comm_cart); MPI_Bcast(&reaction.rate, 1, MPI_DOUBLE, 0, comm_cart); MPI_Bcast(&reaction.back_rate, 1, MPI_DOUBLE, 0, comm_cart); + MPI_Bcast(&reaction.sing_mult, 1, MPI_INT, 0, comm_cart); make_particle_type_exist(reaction.catalyzer_type); make_particle_type_exist(reaction.product_type); @@ -47,7 +48,7 @@ void setup_reaction() { char *error_msg = runtime_error(128); ERROR_SPRINTF(error_msg, "{106 interaction parameters for reaction could not be set} "); } - + data->REACTION_range = reaction.range; /* broadcast interaction parameters */ @@ -55,9 +56,9 @@ void setup_reaction() { } void integrate_reaction() { - int c, np, n, i; - Particle * p1, *p2, **pairs; - Cell * cell; + int c, np, n, i, react; + Particle *p1, *p2, **pairs; + Cell *cell; double dist2, vec21[3], rand; if(reaction.rate > 0) { @@ -85,15 +86,47 @@ void integrate_reaction() { dist2 = sqrlen(vec21); if(dist2 < reaction.range * reaction.range) { - rand = d_random(); - - if(rand > ratexp) { - if(p1->p.type == reaction.reactant_type) { - p1->p.type = reaction.product_type; + + /* Each reactant can react with multiple (neighbouring) catalysts */ + if( reaction.sing_mult == 0 ) + { + rand = d_random(); + + if(rand > ratexp) { + if(p1->p.type == reaction.reactant_type) { + p1->p.type = reaction.product_type; + } + else { + p2->p.type = reaction.product_type; + } + } + } + else /* We only consider each reactant once */ + { + if(p1->p.type == reaction.reactant_type) { + react = p1->p.reacted; } else { - p2->p.type = reaction.product_type; + react = p2->p.reacted; } + + if(react == 0) { + + rand = d_random(); + + if(p1->p.type == reaction.reactant_type) { + if(rand > ratexp) { + p1->p.type = reaction.product_type; + p1->p.reacted = 1; + } + } + else { + if(rand > ratexp) { + p2->p.type = reaction.product_type; + p2->p.reacted = 1; + } + } + } } } } @@ -101,6 +134,19 @@ void integrate_reaction() { } } + /* Clean up the reactants for the next time step */ + if( reaction.sing_mult != 0 ) { + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + p1 = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + if(p1[i].p.type == reaction.product_type) p1[i].p.reacted = 0; + } + } + } + if (reaction.back_rate < 0) { // we have to determine it dynamically /* we count now how many reactants and products are in the sim box */ for (c = 0; c < local_cells.n; c++) { diff --git a/src/reaction.h b/src/reaction.h index 0fb5a865651..c94d80f50b0 100644 --- a/src/reaction.h +++ b/src/reaction.h @@ -34,6 +34,7 @@ typedef struct { double range; double rate; double back_rate; + int sing_mult; } reaction_struct; reaction_struct reaction; diff --git a/src/tcl/initialize_interpreter.c b/src/tcl/initialize_interpreter.c index f79a293c553..3fd45f10d81 100644 --- a/src/tcl/initialize_interpreter.c +++ b/src/tcl/initialize_interpreter.c @@ -28,6 +28,7 @@ #include "constraint_tcl.h" #include "domain_decomposition_tcl.h" #include "dpd_tcl.h" +#include "galilei_tcl.h" #include "global_tcl.h" #include "grid_tcl.h" #include "iccp3m_tcl.h" @@ -207,6 +208,13 @@ static void register_tcl_commands(Tcl_Interp* interp) { #ifdef REACTIONS REGISTER_COMMAND("reaction", tclcommand_reaction); #endif +#ifdef GALILEI + REGISTER_COMMAND("kill_particle_motion", tclcommand_kill_particle_motion); + REGISTER_COMMAND("kill_particle_forces", tclcommand_kill_particle_forces); + REGISTER_COMMAND("system_CMS", tclcommand_system_CMS); + REGISTER_COMMAND("system_CMS_velocity", tclcommand_system_CMS_velocity); + REGISTER_COMMAND("galilei_transform", tclcommand_galilei_transform); +#endif } static void register_global_variables(Tcl_Interp *interp) diff --git a/src/tcl/reaction_tcl.c b/src/tcl/reaction_tcl.c index 4d88e2808db..3d5e672756b 100644 --- a/src/tcl/reaction_tcl.c +++ b/src/tcl/reaction_tcl.c @@ -32,7 +32,7 @@ #ifdef REACTIONS int tcl_command_reaction_print_usage(Tcl_Interp * interp){ char buffer[256]; - sprintf(buffer, "Usage: reaction [off | reactant_type catalyzer_type product_type range rate [back_rate
]]\n"); + sprintf(buffer, "Usage: reaction [off | print | reactant_type catalyzer_type product_type range rate [back_rate
] [react_once ]\n"); Tcl_AppendResult(interp, buffer, (char *)NULL); return TCL_ERROR; } @@ -40,20 +40,25 @@ int tcl_command_reaction_print_usage(Tcl_Interp * interp){ int tcl_command_reaction_print(Tcl_Interp * interp){ char buffer[512]; sprintf(buffer, "{reactant_type %d} {catalyzer_type %d} {product_type %d} {range %f} {rate %f} {back_rate %f}\n", - reaction.reactant_type, - reaction.catalyzer_type, - reaction.product_type, - reaction.range, - reaction.rate, - reaction.back_rate); + reaction.reactant_type, + reaction.catalyzer_type, + reaction.product_type, + reaction.range, + reaction.rate, + reaction.back_rate); + if ( reaction.sing_mult == 0 ) { + sprintf(buffer + strlen(buffer)," {react_once off}\n"); + } else { + sprintf(buffer + strlen(buffer)," {react_once on}\n"); + } Tcl_AppendResult(interp, buffer, (char *)NULL); return TCL_OK; } -#endif /* ifdef REACTIONS */ +#endif -#define ARG_IS_S_EXACT(no, str) !strcmp(argv[(no)], (str)) int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ #ifdef REACTIONS + if (argc == 1 ) return tcl_command_reaction_print_usage(interp); if (argc == 2 ) { if (ARG1_IS_S("off")) { @@ -65,8 +70,9 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** return tcl_command_reaction_print(interp); } } - if( argc!=11 && argc!=13) + if( argc!=11 && argc!=13 && argc!=15) { return tcl_command_reaction_print_usage(interp); + } if(reaction.rate != 0.0) { Tcl_AppendResult(interp, "Currently a simulation can only contain a single reaction!", (char *) NULL); @@ -81,49 +87,57 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** argc--; argv++; while (argc>0){ - if (ARG_IS_S(0,"product_type")) { - if (!ARG_IS_I(1,reaction.product_type)) - return tcl_command_reaction_print_usage(interp); - argc-=2; - argv+=2; - } else - if (ARG_IS_S(0,"reactant_type")) { - if (!ARG_IS_I(1,reaction.reactant_type)) - return tcl_command_reaction_print_usage(interp); - argc-=2; - argv+=2; - } else - if (ARG_IS_S(0,"catalyzer_type")) { - if (!ARG_IS_I(1,reaction.catalyzer_type)) - return tcl_command_reaction_print_usage(interp); - argc-=2; + if (ARG_IS_S(0,"product_type")) { + if (!ARG_IS_I(1,reaction.product_type)) + return tcl_command_reaction_print_usage(interp); + argc-=2; + argv+=2; + } else + if (ARG_IS_S(0,"reactant_type")) { + if (!ARG_IS_I(1,reaction.reactant_type)) + return tcl_command_reaction_print_usage(interp); + argc-=2; + argv+=2; + } else + if (ARG_IS_S(0,"catalyzer_type")) { + if (!ARG_IS_I(1,reaction.catalyzer_type)) + return tcl_command_reaction_print_usage(interp); + argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"range")) { - if (!ARG_IS_D(1,reaction.range)) - return tcl_command_reaction_print_usage(interp); - argc-=2; + } else + if (ARG_IS_S_EXACT(0,"range")) { + if (!ARG_IS_D(1,reaction.range)) + return tcl_command_reaction_print_usage(interp); + argc-=2; + argv+=2; + } else + if (ARG_IS_S_EXACT(0,"rate")) { + if (!ARG_IS_D(1,reaction.rate)) + return tcl_command_reaction_print_usage(interp); + argc-=2; + argv+=2; + } else + if (ARG_IS_S_EXACT(0,"back_rate")) { + if (!ARG_IS_D(1,reaction.back_rate)) + return tcl_command_reaction_print_usage(interp); + argc-=2; + argv+=2; + } else + if (ARG_IS_S_EXACT(0,"react_once")) { + if (!ARG_IS_S(1,"on")&&!ARG_IS_S(1,"off")) { + return tcl_command_reaction_print_usage(interp);} + if (ARG_IS_S(1,"on")) reaction.sing_mult = 1; + if (ARG_IS_S(1,"off")) reaction.sing_mult = 0; + argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"rate")) { - if (!ARG_IS_D(1,reaction.rate)) - return tcl_command_reaction_print_usage(interp); - argc-=2; - argv+=2; - } else - if (ARG_IS_S_EXACT(0,"back_rate")) { - if (!ARG_IS_D(1,reaction.back_rate)) - return tcl_command_reaction_print_usage(interp); - argc-=2; - argv+=2; - } else { - return tcl_command_reaction_print_usage(interp); - } - } - mpi_bcast_event(REACTION); - return TCL_OK; -#else /* ifdef REACTIONS */ + } else { + return tcl_command_reaction_print_usage(interp); + } + } + mpi_bcast_event(REACTION); + return TCL_OK; +#else Tcl_AppendResult(interp, "REACTIONS not compiled in!" ,(char *) NULL); return (TCL_ERROR); -#endif /* ifdef REACTIONS */ +#endif } From ec8db01c54a5069c8bd87f30ef7ece0324165b51 Mon Sep 17 00:00:00 2001 From: "J. de Graaf" Date: Thu, 6 Dec 2012 18:02:45 +0100 Subject: [PATCH 2/5] galilei is all finished up, must still work on the reactions --- scripts/auxiliary.tcl | 208 ++++++++++++++++++++--------------------- src/communication.c | 140 +++++++++++++++++++++++++++ src/communication.h | 8 ++ src/tcl/reaction_tcl.c | 26 +++--- 4 files changed, 263 insertions(+), 119 deletions(-) diff --git a/scripts/auxiliary.tcl b/scripts/auxiliary.tcl index e4a5ade9bdc..68e06d5a0f8 100644 --- a/scripts/auxiliary.tcl +++ b/scripts/auxiliary.tcl @@ -445,24 +445,24 @@ proc tune_cells { { int_steps 1000 } { min_cells 1 } { max_cells 30 } { tol_cell # stored in Espresso to zero. # ############################################################# - -proc stop_particles { } { - for { set i 0} { $i < [setmd n_part] } {incr i} { - part $i v 0 0 0 - part $i f 0 0 0 - } -} -proc stopParticles { } { - puts -nonewline " Setting all particles' velocities and forces to zero... "; flush stdout - set old_i 0; set old_e [setmd n_part]; set old [expr 10*$old_e] - for { set i 0} { $i < $old_e } {incr i} { - if { [expr $old_i*100]>=$old } { set old [expr $old + 10*$old_e]; puts -nonewline ".,."; flush stdout }; incr old_i - part $i v 0 0 0 - part $i f 0 0 0 - } - puts ".,. Done (stopped $old_e particles and as many forces)." -} - +# +#proc stop_particles { } { +# for { set i 0} { $i < [setmd n_part] } {incr i} { +# part $i v 0 0 0 +# part $i f 0 0 0 +# } +#} +#proc stopParticles { } { +# puts -nonewline " Setting all particles' velocities and forces to zero... "; flush stdout +# set old_i 0; set old_e [setmd n_part]; set old [expr 10*$old_e] +# for { set i 0} { $i < $old_e } {incr i} { +# if { [expr $old_i*100]>=$old } { set old [expr $old + 10*$old_e]; puts -nonewline ".,."; flush stdout }; incr old_i +# part $i v 0 0 0 +# part $i f 0 0 0 +# } +# puts ".,. Done (stopped $old_e particles and as many forces)." +#} +# # # center of mass # -------------- @@ -470,41 +470,37 @@ proc stopParticles { } { # Calculate the center of mass of the system stored in Espresso # ############################################################# -proc system_com { } { - set npart [setmd n_part] - - #calculate center of mass - set com {0.0 0.0 0.0} - set part_cnt 0 - - if {[has_feature MASS]} { - set net_mass 0 - - for {set i 0} {$i < $npart} {incr i} { - if {[part $i] != "na"} { - set pos [part $i print p] - set mass [part $i print mass] - set com [vecadd $com [vecscale $mass $pos]] - set net_mass [expr $net_mass + $mass] - incr part_cnt - } - } - - return [vecscale [expr 1.0/$net_mass] $com] - } else { - for {set i 0} {$i < $npart} {incr i} { - if {[part $i] != "na"} { - set pos [part $i print p] - set com [vecadd $com $pos] - incr part_cnt - } - } - - return [vecscale [expr 1.0/$npart] $com] - } -} - - +#proc system_com { } { +# set npart [setmd n_part] +# +# #calculate center of mass +# set com {0.0 0.0 0.0} +# set part_cnt 0 +# +# if {[has_feature MASS]} { +# set net_mass 0 +# +# for {set i 0} {$i < $npart} {incr i} { +# set pos [part $i print p] +# set mass [part $i print mass] +# set com [vecadd $com [vecscale $mass $pos]] +# set net_mass [expr $net_mass + $mass] +# incr part_cnt +# } +# +# return [vecscale [expr 1.0/$net_mass] $com] +# } else { +# for {set i 0} {$i < $npart} {incr i} { +# set pos [part $i print p] +# set com [vecadd $com $pos] +# incr part_cnt +# } +# +# return [vecscale [expr 1.0/$npart] $com] +# } +#} +# +# # # center of mass motion # --------------------- @@ -512,40 +508,40 @@ proc system_com { } { # Calculate the center of mass velocity of the system stored in Espresso # ############################################################# -proc system_com_vel {} { - set npart [setmd n_part] - - # calculate center of mass motion - set com_vel {0 0 0} - set part_cnt 0 - - if {[has_feature MASS]} { - set net_mass 0 - - for {set i 0} {$i < $npart} {incr i} { - if { [part $i] != "na" } { - set part_vel [part $i print v] - set part_mass [part $i print mass] - set com_vel [vecadd $com_vel [vecscale $part_mass $part_vel]] - set net_mass [expr $net_mass + $part_mass] - incr part_cnt - } - } - - return "[vecscale [expr 1.0/$net_mass] $com_vel]" - } else { - for {set i 0} {$i < $npart} {incr i} { - if { [part $i] != "na" } { - set part_vel [part $i print v] - set com_vel [vecadd $com_vel $part_vel] - incr part_cnt - } - } - - return [vecscale [expr 1.0/$npart] $com_vel] - } -} - +#proc system_com_vel {} { +# set npart [setmd n_part] +# +# # calculate center of mass motion +# set com_vel {0 0 0} +# set part_cnt 0 +# +# if {[has_feature MASS]} { +# set net_mass 0 +# +# for {set i 0} {$i < $npart} {incr i} { +# if { [part $i] != "na" } { +# set part_vel [part $i print v] +# set part_mass [part $i print mass] +# set com_vel [vecadd $com_vel [vecscale $part_mass $part_vel]] +# set net_mass [expr $net_mass + $part_mass] +# incr part_cnt +# } +# } +# +# return "[vecscale [expr 1.0/$net_mass] $com_vel]" +# } else { +# for {set i 0} {$i < $npart} {incr i} { +# if { [part $i] != "na" } { +# set part_vel [part $i print v] +# set com_vel [vecadd $com_vel $part_vel] +# incr part_cnt +# } +# } +# +# return [vecscale [expr 1.0/$npart] $com_vel] +# } +#} +# # # galileiTransformParticles # ------------------------- @@ -555,25 +551,25 @@ proc system_com_vel {} { # system is zero. Returns old center of mass motion. # ############################################################# - -proc galileiTransformParticles {} { - set com_vel [system_com_vel] - - #subtract center of mass motion from particle velocities - set npart [setmd n_part]; - set part_cnt 0; - - for {set i 0} { $part_cnt < $npart } {incr i} { - if {[part $i] != "na"} { - set part_vel [vecsub [part $i print v] $com_vel] - part $i v [lindex $part_vel 0] [lindex $part_vel 1] [lindex $part_vel 2] - incr part_cnt - } - } - - return $com_vel -} - +# +#proc galileiTransformParticles {} { +# set com_vel [system_com_vel] +# +# #subtract center of mass motion from particle velocities +# set npart [setmd n_part]; +# set part_cnt 0; +# +# for {set i 0} { $part_cnt < $npart } {incr i} { +# if {[part $i] != "na"} { +# set part_vel [vecsub [part $i print v] $com_vel] +# part $i v [lindex $part_vel 0] [lindex $part_vel 1] [lindex $part_vel 2] +# incr part_cnt +# } +# } +# +# return $com_vel +#} +# # # Prepare Connection to VMD # ------------------------- diff --git a/src/communication.c b/src/communication.c index 622af62d46e..25216b581d7 100644 --- a/src/communication.c +++ b/src/communication.c @@ -61,6 +61,7 @@ #include "molforces.h" #include "mdlc_correction.h" #include "reaction.h" +#include "galilei.h" int this_node = -1; int n_nodes = -1; @@ -137,6 +138,11 @@ typedef void (SlaveCallback)(int node, int param); CB(mpi_recv_fluid_boundary_flag_slave) \ CB(mpi_set_particle_temperature_slave) \ CB(mpi_set_particle_gamma_slave) \ + CB(mpi_kill_particle_motion_slave) \ + CB(mpi_kill_particle_forces_slave) \ + CB(mpi_system_CMS_slave) \ + CB(mpi_system_CMS_velocity_slave) \ + CB(mpi_galilei_transform_slave) \ // create the forward declarations #define CB(name) void name(int node, int param); @@ -2637,6 +2643,140 @@ void mpi_set_particle_gamma_slave(int pnode, int part) #endif } +/******************** REQ_GALILEI ********************/ + +void mpi_kill_particle_motion() +{ + mpi_call(mpi_kill_particle_motion_slave, -1, 0); + local_kill_particle_motion(); + on_particle_change(); +} + +void mpi_kill_particle_motion_slave(int pnode, int i) +{ + local_kill_particle_motion(); + on_particle_change(); +} + +void mpi_kill_particle_forces() +{ + mpi_call(mpi_kill_particle_forces_slave, -1, 0); + local_kill_particle_forces(); + on_particle_change(); +} + +void mpi_kill_particle_forces_slave(int pnode, int i) +{ + local_kill_particle_forces(); + on_particle_change(); +} + +void mpi_system_CMS() { + int pnode; + double data[4]; + double rdata[4]; + double *pdata = rdata; + + data[0] = 0.0; + data[1] = 0.0; + data[2] = 0.0; + data[3] = 0.0; + + mpi_call(mpi_system_CMS_slave, -1, 0); + + for (pnode = 0; pnode < n_nodes; pnode++) { + if (pnode==this_node) { + local_system_CMS( pdata ); + data[0] += rdata[0]; + data[1] += rdata[1]; + data[2] += rdata[2]; + data[3] += rdata[3]; + } else { + MPI_Recv(rdata, 4, MPI_DOUBLE, MPI_ANY_SOURCE, SOME_TAG, comm_cart, MPI_STATUS_IGNORE); + data[0] += rdata[0]; + data[1] += rdata[1]; + data[2] += rdata[2]; + data[3] += rdata[3]; + } + } + + gal.cms[0] = data[0]/data[3]; + gal.cms[1] = data[1]/data[3]; + gal.cms[2] = data[2]/data[3]; +} + +void mpi_system_CMS_slave(int node, int index) { + double rdata[4]; + double *pdata = rdata; + local_system_CMS( pdata ); + MPI_Send(rdata, 4, MPI_DOUBLE, 0, SOME_TAG, comm_cart); +} + +void mpi_system_CMS_velocity() { + int pnode; + double data[4]; + double rdata[4]; + double *pdata = rdata; + + data[0] = 0.0; + data[1] = 0.0; + data[2] = 0.0; + data[3] = 0.0; + + mpi_call(mpi_system_CMS_velocity_slave, -1, 0); + + for (pnode = 0; pnode < n_nodes; pnode++) { + if (pnode==this_node) { + local_system_CMS_velocity( pdata ); + data[0] += rdata[0]; + data[1] += rdata[1]; + data[2] += rdata[2]; + data[3] += rdata[3]; + } else { + MPI_Recv(rdata, 4, MPI_DOUBLE, MPI_ANY_SOURCE, SOME_TAG, comm_cart, MPI_STATUS_IGNORE); + data[0] += rdata[0]; + data[1] += rdata[1]; + data[2] += rdata[2]; + data[3] += rdata[3]; + } + } + + gal.cms_vel[0] = data[0]/data[3]; + gal.cms_vel[1] = data[1]/data[3]; + gal.cms_vel[2] = data[2]/data[3]; +} + +void mpi_system_CMS_velocity_slave(int node, int index) { + double rdata[4]; + double *pdata = rdata; + local_system_CMS_velocity( pdata ); + MPI_Send(rdata, 4, MPI_DOUBLE, 0, SOME_TAG, comm_cart); +} + +void mpi_galilei_transform() +{ + double cmsvel[3]; + + mpi_system_CMS_velocity(); + memcpy(cmsvel, gal.cms_vel, 3*sizeof(double)); + + mpi_call(mpi_galilei_transform_slave, -1, 0); + MPI_Bcast(cmsvel, 3, MPI_DOUBLE, 0, comm_cart); + + local_galilei_transform( cmsvel ); + + on_particle_change(); +} + +void mpi_galilei_transform_slave(int pnode, int i) +{ + double cmsvel[3]; + MPI_Bcast(cmsvel, 3, MPI_DOUBLE, 0, comm_cart); + + local_galilei_transform( cmsvel ); + on_particle_change(); +} + /*********************** MAIN LOOP for slaves ****************/ void mpi_loop() diff --git a/src/communication.h b/src/communication.h index 4736d4a4dd5..3759f6fb21c 100644 --- a/src/communication.h +++ b/src/communication.h @@ -518,6 +518,14 @@ void mpi_bcast_max_mu(); */ int mpi_gather_runtime_errors(char **errors); +/** Issue REQ_GALILEI: set all particle velocities and rotational inertias to zero. + set all forces and torques on the particles to zero */ +void mpi_kill_particle_motion(); +void mpi_kill_particle_forces(); +void mpi_system_CMS(); +void mpi_system_CMS_velocity(); +void mpi_galilei_transform(); + /*@}*/ /** \name Event codes for \ref mpi_bcast_event diff --git a/src/tcl/reaction_tcl.c b/src/tcl/reaction_tcl.c index 3d5e672756b..3c4c9f9cf48 100644 --- a/src/tcl/reaction_tcl.c +++ b/src/tcl/reaction_tcl.c @@ -70,7 +70,7 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** return tcl_command_reaction_print(interp); } } - if( argc!=11 && argc!=13 && argc!=15) { + if( argc!=11 && argc!=13 && argc!=15 && argc!=3) { return tcl_command_reaction_print_usage(interp); } @@ -92,38 +92,38 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S(0,"reactant_type")) { + } + else if (ARG_IS_S(0,"reactant_type")) { if (!ARG_IS_I(1,reaction.reactant_type)) return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S(0,"catalyzer_type")) { + } + else if (ARG_IS_S(0,"catalyzer_type")) { if (!ARG_IS_I(1,reaction.catalyzer_type)) return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"range")) { + } + else if (ARG_IS_S_EXACT(0,"range")) { if (!ARG_IS_D(1,reaction.range)) return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"rate")) { + } + else if (ARG_IS_S_EXACT(0,"rate")) { if (!ARG_IS_D(1,reaction.rate)) return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"back_rate")) { + } + else if (ARG_IS_S_EXACT(0,"back_rate")) { if (!ARG_IS_D(1,reaction.back_rate)) return tcl_command_reaction_print_usage(interp); argc-=2; argv+=2; - } else - if (ARG_IS_S_EXACT(0,"react_once")) { + } + else if (ARG_IS_S_EXACT(0,"react_once")) { if (!ARG_IS_S(1,"on")&&!ARG_IS_S(1,"off")) { return tcl_command_reaction_print_usage(interp);} if (ARG_IS_S(1,"on")) reaction.sing_mult = 1; From c19553c37edf39078faa9c294509353657a47604 Mon Sep 17 00:00:00 2001 From: "J. de Graaf" Date: Tue, 8 Jan 2013 11:43:34 +0100 Subject: [PATCH 3/5] This commit updates the following elements of ESPResSo. 1) It removes the Tcl functions: stop_particles, stopParticles, system_com, system_com_vel, and galileiTransformParticles, which were in the auxiliary.tcl script. 2) It replaces these with MPI functions: kill_particle_motion, kill_particle_forces, system_CMS, system_CMS_velocity, and galilei_transform (respectively). These new functions are incorporated in a new feature called 'GALILEI'. 3) It changes the Hat-potential entry in the user guide to be more clear on what the cut-off radius is and what the radial dependence is. 4) It updates the reaction code with two options: (a) Letting each reactant interact only with one catalyst particle during each time step. This enables the creation of catalytic surfaces comprised of catalyst particles (with arbitrary density), yet maintaining a surface-density-independent reaction rate near the surface. (b) A print option is given to the reaction command, such that the specific input parameters can be returned to the Tcl level. 5) All the necessary updates are made in the user guide. --- doc/ug/aux.tex | 28 ------------- doc/ug/features.tex | 1 + doc/ug/inter.tex | 23 +++++------ doc/ug/run.tex | 14 +++---- doc/ug/setup.tex | 91 +++++++++++++++++++++++++++++++----------- src/communication.c | 43 +++++++++++++++++--- src/communication.h | 12 +++++- src/features.def | 1 + src/reaction.c | 38 +++++++++--------- src/reaction.h | 2 +- src/tcl/reaction_tcl.c | 4 +- 11 files changed, 156 insertions(+), 101 deletions(-) diff --git a/doc/ug/aux.tex b/doc/ug/aux.tex index 55df0fed716..38f26317188 100644 --- a/doc/ug/aux.tex +++ b/doc/ug/aux.tex @@ -23,34 +23,6 @@ \chapter{Auxilliary commands} \todo{Missing commands: Probably all from \keyword{scripts/auxiliary.tcl}?} -\section{Center of mass motion} -\label{sec:centerofmassmotion} - -\subsection{\texttt{system_com}} -\begin{essyntax} - system_com -\end{essyntax} -Returns the center of mass of the whole system. - -\subsection{\texttt{system_com_vel}} -\begin{essyntax} - system_com_vel -\end{essyntax} -Returns the velocity of the center of mass of the whole system. - -\subsection{\texttt{galileiTransformParticles}} -\begin{essyntax} - galileiTransformParticles -\end{essyntax} -Substracts the velocity of the center of mass of the whole system from every -particle's velocity, thereby performing a galilei transform into the reference -frame of the center of mass of the system. - -This is useful for example in combination with the DPD thermostat, since there, - a drift in the velocity of the whole system leads to an offset in the reported - temperature. - - \section{Finding particles and bonds} \subsection{\texttt{countBonds}} diff --git a/doc/ug/features.tex b/doc/ug/features.tex index e230ceb8bb8..00d6b4aa58b 100644 --- a/doc/ug/features.tex +++ b/doc/ug/features.tex @@ -112,6 +112,7 @@ \section{General features} \item \newfeature{OLD\_RW\_VERSION} This switches back to the old, \emph{wrong} random walk code of the polymer. Only use this if you rely on the old behaviour and \emph{know what you are doing}. +\item \newfeature{GALILEI} Allows for a global Galilei transform of the system. It also holds routines to determine the system's center of mass and it's velocity, as well as routines to stop the particles from moving and set the forces on all particles to zero. \end{itemize} In addition, there are switches that enable additional features in the diff --git a/doc/ug/inter.tex b/doc/ug/inter.tex index 1747501167c..b8d3dffcb8a 100644 --- a/doc/ug/inter.tex +++ b/doc/ug/inter.tex @@ -421,32 +421,31 @@ \subsection{Hat interaction} \begin{essyntax} inter \var{type1} \var{type2} - hat \var{F_\text{max}} \var{r} + hat \var{F_\text{max}} \var{r_c} \begin{features} \required{HAT} \end{features} \end{essyntax} This defines a simple force ramp between particles of the types \var{type1} and \var{type2}. The maximal force \var{F_\text{max}} acts at zero distance and -zero force is applied at distances r and bigger. For distances smaller than -\var{r}, the force is given by +zero force is applied at distances $r_c$ and bigger. For distances smaller than +\var{r_c}, the force is given by \begin{equation} - F(d)=F_{\text{max}} \cdot \left( 1 - \frac{d}{r} \right), + F(r)=F_{\text{max}} \cdot \left( 1 - \frac{r}{r_c} \right), \end{equation} -for distances exceeding \var{r}, the force is zero. +for distances exceeding \var{r_c}, the force is zero. The potential energy is given by \begin{equation} - V(d)=F_{\text{max}} \cdot (d-r) \cdot \left( \frac{d+r}{2r} - 1 \right), + V(r)=F_{\text{max}} \cdot (r-r_c) \cdot \left( \frac{r+r_c}{2r_c} - 1 \right), \end{equation} -which is zero for distances bigger than \var{r} and continuous at distance -\var{r}. +which is zero for distances bigger than \var{r_c} and continuous at distance +\var{r_c}. This is the standard conservative DPD potential and can be used in combination -with \lit{inter DPD} \ref{sec:DPDinter}. Also it is useful for live -demonstrations, where you might want to use a big time step to make it look -interesting on a weak machine but do not want it to crash and the physics don't -matter. +with \lit{inter DPD} \ref{sec:DPDinter}. The potential is also useful for live +demonstrations, where a big time step may be employed to obtain quick results +on a weak machine, for which the physics do not need to be entirely correct. \subsection{Hertzian interaction} \index{Hertzian interaction|mainindex} diff --git a/doc/ug/run.tex b/doc/ug/run.tex index d20da6055d7..2ca8d0b41c4 100644 --- a/doc/ug/run.tex +++ b/doc/ug/run.tex @@ -64,16 +64,12 @@ \section{\texttt{change_volume}: Changing the box volume} the deformed simulation box. \section{Stopping particles} -\newescommand{stopParticles} -\newescommand[stop-particles]{stop_particles} -\begin{essyntax} - \variant{1} stopParticles - \variant{2} stop_particles -\end{essyntax} -Halts all particles in the current simulation, setting their -velocities and forces to zero. Variant \variant{2} does not provide -feedback on the execution status. +Use the following functions that compile in with the GALILEI feature \ref{sec:Galilei}: +\begin{itemize} +\item \texttt{kill\_particle\_motion}: halts all particles in the current simulation, setting their velocities to zero, as well as their angular momentum if the feature ROTATION has been compiled in. +\item \texttt{kill\_particle\_forces}: sets all forces on the particles to zero, as well as all torques if the feature ROTATION has been compiled in. +\end{itemize} \section{\texttt{velocities}: Setting the velocities} \newescommand{velocities} diff --git a/doc/ug/setup.tex b/doc/ug/setup.tex index e9e465d1367..ac4512aa133 100644 --- a/doc/ug/setup.tex +++ b/doc/ug/setup.tex @@ -306,10 +306,10 @@ \subsubsection{Thermostat DPD} temperature jumps during a simulation are. The dpd thermostat does not act on the system center of mass motion. Therefore, before using dpd, you have to stop the center of mass motion of your system, which you -can achieve by using the command \keyword{galileiTransformParticles} -\ref{sec:centerofmassmotion}. This may be repeated once in a while for long +can achieve by using the command \keyword{galilei_transform} +\ref{sec:Galilei}. This may be repeated once in a while for long runs due to round off errors (check this with the command -\keyword{system_com_vel}) \ref{sec:centerofmassmotion}. +\keyword{system_CMS_velocity}) \ref{sec:Galilei}. Two restrictions apply for the dpd implementation of \es: @@ -679,40 +679,85 @@ \section{Creating bonds when particles collide} \section{Reactions} \label{sec:Reactions} -With the help of the feature \feature{REACTIONS}, one can define three particle -types to act as reactant, catalyzer and product. Reactant particles that are -closer to at least one catalyzer particle than a user defined range are -stochastically converted into particles of the product type. The rate at which -this reaction happens can be defined by the user as well. Since an Espresso -simulation only contains a finite number of particles, product particles also -must react back into reactant particles in order allow for an equilibrium other -than only product particles. -This back reaction takes place at any point of the simulation volume and unless -the user defines a custom rate for the back reaction, the rate is dynamically -chosen to be +With the help of the feature \feature{REACTIONS}, one can define three particle types to act as reactant, catalyzer and product. Reactant particles that are closer to at least one catalyzer particle than a user defined range are stochastically converted into particles of the product type. The rate at which this reaction happens can be defined by the user as well. Since an \es simulation only contains a finite number of particles, product particles also +must react back into reactant particles in order allow for an equilibrium other than only product particles. This back reaction takes place at any point of the simulation volume and unless +the user defines a custom rate for the back reaction, the rate is dynamically chosen to be \begin{equation} k_\text{back} = k_\text{back} \cdot \frac{n_\text{product}}{n_\text{reactant}}, \end{equation} -which results in an equilibrium of half product and half reactant, independently -of the number of catalyzer particles or their arrangement. +which results in an equilibrium of half product and half reactant, independently of the number of catalyzer particles or their arrangement. \begin{essyntax} \variant{0} reaction reactant_type \var{rt} catalyzer_type \var{ct} product_type \var{pt} range \var{r} rate \var{k} -\opt{back_rate \var{k_\text{back}}} +\opt{back_rate \var{k_\text{back}}} \opt{react_once \var{on/off}} \variant{1} reaction off +\variant{2} reaction print \begin{features} \required{REACTIONS} \end{features} \end{essyntax} -Variant \variant{0} defines a reaction with particles of type number \var{rt} -as reactant, type \var{ct} as catalyzer and type \var{pt} as product. The -reaction rate is given by \var{r} and to override the dynamically determined -rate for the back reaction, one can specify it by \var{k_\text{back}}. +\begin{itemize} +\item Variant \variant{0} defines a reaction with particles of type number \var{rt} as reactant, type \var{ct} as catalyzer and type \var{pt} as product. The reaction rate is given by \var{r} and to override the dynamically determined rate for the back reaction, one can specify it by \var{k_\text{back}}. By default each reactant particle is checked against each catalyst particle (\texttt{react\_once} \emph{off}). However, when creating smooth surfaces using many catalyst particles, it can be desirable to let the reaction rate be independent of the surface density of these particles. That is, each particle has a likelyhood of reacting in the vacinity of the surface (distance is less than $r$) as specified by the rate, \emph{not} by the rate muliplied by the local catalyst surface density. To accomplish this, each reactant is considered only once each time step by using the option \texttt{react\_once} \emph{on}. The reaction command is set up such that the different properties may be influenced individually. +\item Variant \variant{1} disables the reaction. Note that at the moment, there can only be one reaction in the simulation. +\item Variant \variant{2} returns the current reaction parameters. +\end{itemize} + +\section{Galilei} +\label{sec:Galilei} + +The GALILEI feature holds several commands which may be useful in effecting the velocity of the system. + +\subsection{Particle motion and rotation} + +\begin{essyntax} + kill_particle_motion +\begin{features} + \required{GALILEI} +\end{features} +\end{essyntax} +This command halts all particles in the current simulation, setting their velocities to zero, as well as their angular momentum if the feature ROTATION has been compiled in. + +\subsection{Forces and torques acting on the particles} + +\begin{essyntax} + kill_particle_forces +\begin{features} + \required{GALILEI} +\end{features} +\end{essyntax} +This command sets all forces on the particles to zero, as well as all torques if the feature ROTATION has been compiled in. + +\subsection{The centre of mass of the system} -Variant \variant{1} disables the reaction. Note that at the moment, there can -only be one reaction in the simulation. +\begin{essyntax} + system_CMS +\begin{features} + \required{GALILEI} +\end{features} +\end{essyntax} +Returns the center of mass of the whole system. It currently does not factor in the density fluctuations of the Lattice-Boltzman fluid. + +\subsection{The centre-of-mass velocity} + +\begin{essyntax} + system_CMS_velocity +\begin{features} + \required{GALILEI} +\end{features} +\end{essyntax} +Returns the velocity of the center of mass of the whole system. + +\subsection{The Galilei transform} + +\begin{essyntax} + galilei_transform +\begin{features} + \required{GALILEI} +\end{features} +\end{essyntax} +Substracts the velocity of the center of mass of the whole system from every particle's velocity, thereby performing a Galilei transform into the reference frame of the center of mass of the system. This transformation is useful for example in combination with the DPD thermostat, since there, a drift in the velocity of the whole system leads to an offset in the reported temperature. %%% Local Variables: %%% mode: latex diff --git a/src/communication.c b/src/communication.c index 25216b581d7..07f8ce515b7 100644 --- a/src/communication.c +++ b/src/communication.c @@ -143,6 +143,7 @@ typedef void (SlaveCallback)(int node, int param); CB(mpi_system_CMS_slave) \ CB(mpi_system_CMS_velocity_slave) \ CB(mpi_galilei_transform_slave) \ + CB(mpi_setup_reaction_slave) \ // create the forward declarations #define CB(name) void name(int node, int param); @@ -405,11 +406,6 @@ void mpi_bcast_event_slave(int node, int event) break; #endif -#ifdef REACTIONS - case REACTION: - setup_reaction(); -#endif - default:; } } @@ -2647,31 +2643,40 @@ void mpi_set_particle_gamma_slave(int pnode, int part) void mpi_kill_particle_motion() { +#ifdef GALILEI mpi_call(mpi_kill_particle_motion_slave, -1, 0); local_kill_particle_motion(); on_particle_change(); +#endif } void mpi_kill_particle_motion_slave(int pnode, int i) { +#ifdef GALILEI local_kill_particle_motion(); on_particle_change(); +#endif } void mpi_kill_particle_forces() { +#ifdef GALILEI mpi_call(mpi_kill_particle_forces_slave, -1, 0); local_kill_particle_forces(); on_particle_change(); +#endif } void mpi_kill_particle_forces_slave(int pnode, int i) { +#ifdef GALILEI local_kill_particle_forces(); on_particle_change(); +#endif } void mpi_system_CMS() { +#ifdef GALILEI int pnode; double data[4]; double rdata[4]; @@ -2703,16 +2708,20 @@ void mpi_system_CMS() { gal.cms[0] = data[0]/data[3]; gal.cms[1] = data[1]/data[3]; gal.cms[2] = data[2]/data[3]; +#endif } void mpi_system_CMS_slave(int node, int index) { +#ifdef GALILEI double rdata[4]; double *pdata = rdata; local_system_CMS( pdata ); MPI_Send(rdata, 4, MPI_DOUBLE, 0, SOME_TAG, comm_cart); +#endif } void mpi_system_CMS_velocity() { +#ifdef GALILEI int pnode; double data[4]; double rdata[4]; @@ -2744,17 +2753,21 @@ void mpi_system_CMS_velocity() { gal.cms_vel[0] = data[0]/data[3]; gal.cms_vel[1] = data[1]/data[3]; gal.cms_vel[2] = data[2]/data[3]; +#endif } void mpi_system_CMS_velocity_slave(int node, int index) { +#ifdef GALILEI double rdata[4]; double *pdata = rdata; local_system_CMS_velocity( pdata ); MPI_Send(rdata, 4, MPI_DOUBLE, 0, SOME_TAG, comm_cart); +#endif } void mpi_galilei_transform() { +#ifdef GALILEI double cmsvel[3]; mpi_system_CMS_velocity(); @@ -2766,15 +2779,35 @@ void mpi_galilei_transform() local_galilei_transform( cmsvel ); on_particle_change(); +#endif } void mpi_galilei_transform_slave(int pnode, int i) { +#ifdef GALILEI double cmsvel[3]; MPI_Bcast(cmsvel, 3, MPI_DOUBLE, 0, comm_cart); local_galilei_transform( cmsvel ); on_particle_change(); +#endif +} + +/******************** REQ_REACTIONS ********************/ + +void mpi_setup_reaction() +{ +#ifdef REACTIONS + mpi_call(mpi_setup_reaction_slave, -1, 0); + local_setup_reaction(); +#endif +} + +void mpi_setup_reaction_slave(int pnode, int i) +{ +#ifdef REACTIONS + local_setup_reaction(); +#endif } /*********************** MAIN LOOP for slaves ****************/ diff --git a/src/communication.h b/src/communication.h index 3759f6fb21c..93e3f26d9ec 100644 --- a/src/communication.h +++ b/src/communication.h @@ -519,13 +519,21 @@ void mpi_bcast_max_mu(); int mpi_gather_runtime_errors(char **errors); /** Issue REQ_GALILEI: set all particle velocities and rotational inertias to zero. - set all forces and torques on the particles to zero */ + set all forces and torques on the particles to zero + calculate the centre of mass (CMS) + calculate the velocity of the CMS + remove the CMS velocity from the system + */ void mpi_kill_particle_motion(); void mpi_kill_particle_forces(); void mpi_system_CMS(); void mpi_system_CMS_velocity(); void mpi_galilei_transform(); +/** Issue REQ_REACTIONS: notify the system of changes to the reaction parameters + */ +void mpi_setup_reaction(); + /*@}*/ /** \name Event codes for \ref mpi_bcast_event @@ -538,7 +546,7 @@ void mpi_galilei_transform(); #define CHECK_PARTICLES 2 #define MAGGS_COUNT_CHARGES 3 #define P3M_COUNT_DIPOLES 5 -#define REACTION 6 +//#define REACTION 6 /*@}*/ #endif diff --git a/src/features.def b/src/features.def index a2843b407be..54ceaf1fa7b 100644 --- a/src/features.def +++ b/src/features.def @@ -27,6 +27,7 @@ NEMD NPT GHMC REACTIONS +GALILEI /* Rotation */ ROTATION diff --git a/src/reaction.c b/src/reaction.c index 2c800bc6b16..8bf5376674d 100644 --- a/src/reaction.c +++ b/src/reaction.c @@ -29,7 +29,7 @@ #ifdef REACTIONS -void setup_reaction() { +void local_setup_reaction() { MPI_Bcast(&reaction.reactant_type, 1, MPI_INT, 0, comm_cart); MPI_Bcast(&reaction.product_type, 1, MPI_INT, 0, comm_cart); MPI_Bcast(&reaction.catalyzer_type, 1, MPI_INT, 0, comm_cart); @@ -105,28 +105,28 @@ void integrate_reaction() { { if(p1->p.type == reaction.reactant_type) { react = p1->p.reacted; + p1->p.reacted = 1; + + if(react == 0) { + + rand = d_random(); + if(rand > ratexp) { + p1->p.type = reaction.product_type; + } + } } else { react = p2->p.reacted; - } + p2->p.reacted = 1; - if(react == 0) { - - rand = d_random(); - - if(p1->p.type == reaction.reactant_type) { - if(rand > ratexp) { - p1->p.type = reaction.product_type; - p1->p.reacted = 1; - } - } - else { - if(rand > ratexp) { + if(react == 0) { + + rand = d_random(); + if(rand > ratexp) { p2->p.type = reaction.product_type; - p2->p.reacted = 1; } - } - } + } + } } } } @@ -142,7 +142,7 @@ void integrate_reaction() { np = cell->n; for(i = 0; i < np; i++) { - if(p1[i].p.type == reaction.product_type) p1[i].p.reacted = 0; + if(p1[i].p.type == reaction.reactant_type || p1[i].p.type == reaction.product_type ) p1[i].p.reacted = 0; } } } @@ -192,4 +192,4 @@ void integrate_reaction() { on_particle_change(); } } -#endif /* ifdef REACTIONS */ +#endif diff --git a/src/reaction.h b/src/reaction.h index c94d80f50b0..391b5254630 100644 --- a/src/reaction.h +++ b/src/reaction.h @@ -43,7 +43,7 @@ reaction_struct reaction; /** broadcasts reaction parameters and sets up an entry in the ia_params, so that the verlet radius is equal or bigger than the reaction range. **/ -void setup_reaction(); +void local_setup_reaction(); void integrate_reaction(); #endif diff --git a/src/tcl/reaction_tcl.c b/src/tcl/reaction_tcl.c index 3c4c9f9cf48..4495d89a42f 100644 --- a/src/tcl/reaction_tcl.c +++ b/src/tcl/reaction_tcl.c @@ -63,7 +63,7 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** if (argc == 2 ) { if (ARG1_IS_S("off")) { reaction.rate=0.0; - mpi_bcast_event(REACTION); + mpi_setup_reaction(); return TCL_OK; } if (ARG1_IS_S("print")) { @@ -134,7 +134,7 @@ int tclcommand_reaction(ClientData data, Tcl_Interp * interp, int argc, char ** return tcl_command_reaction_print_usage(interp); } } - mpi_bcast_event(REACTION); + mpi_setup_reaction(); return TCL_OK; #else Tcl_AppendResult(interp, "REACTIONS not compiled in!" ,(char *) NULL); From 83052796e1c6a9b1adca1986ec245a8ad327b970 Mon Sep 17 00:00:00 2001 From: "J. de Graaf" Date: Tue, 8 Jan 2013 12:05:40 +0100 Subject: [PATCH 4/5] Merge branches 'galilei' and 'galilei_transform' into galilei_transform From 540498960c712111a50c3cd57700ead99d54fc23 Mon Sep 17 00:00:00 2001 From: "J. de Graaf" Date: Tue, 8 Jan 2013 12:15:04 +0100 Subject: [PATCH 5/5] Also included the files I missed earlier. --- src/galilei.c | 235 ++++++++++++++++++++++++++++++++++++++++++ src/galilei.h | 52 ++++++++++ src/tcl/galilei_tcl.c | 216 ++++++++++++++++++++++++++++++++++++++ src/tcl/galilei_tcl.h | 45 ++++++++ 4 files changed, 548 insertions(+) create mode 100644 src/galilei.c create mode 100644 src/galilei.h create mode 100644 src/tcl/galilei_tcl.c create mode 100644 src/tcl/galilei_tcl.h diff --git a/src/galilei.c b/src/galilei.c new file mode 100644 index 00000000000..cebab290d7a --- /dev/null +++ b/src/galilei.c @@ -0,0 +1,235 @@ +/* + Copyright (C) 2010,2011,2012 The ESPResSo project + Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + Max-Planck-Institute for Polymer Research, Theory Group + + This file is part of ESPResSo. + + ESPResSo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ESPResSo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +/** \file galilei.c + * + */ + +#include "galilei.h" +#include "utils.h" +#include "initialize.h" +#include "forces.h" + +#ifdef GALILEI + +galilei_struct gal; + +/* Stop the particle motion by setting the + velocity of each particle to zero */ +void local_kill_particle_motion() { + int c, np, i; + Particle *part; + Cell *cell; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + part[i].m.v[0] = 0.0; + part[i].m.v[1] = 0.0; + part[i].m.v[2] = 0.0; + +#ifdef ROTATION + part[i].m.omega[0] = 0.0; + part[i].m.omega[1] = 0.0; + part[i].m.omega[2] = 0.0; +#endif + } + } +} + +/* Set all the forces acting on the particles + to zero */ +void local_kill_particle_forces() { + int c, np, i; + Particle *part; + Cell *cell; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + part[i].f.f[0] = 0.0; + part[i].f.f[1] = 0.0; + part[i].f.f[2] = 0.0; + +#ifdef ROTATION + part[i].f.torque[0] = 0.0; + part[i].f.torque[1] = 0.0; + part[i].f.torque[2] = 0.0; +#endif + } + } +} + +/* Calculate the CMS of the system */ +void local_system_CMS( double *sdata ) { + int c, np, i; + Particle *part; + Cell *cell; + double x = 0.0, + y = 0.0, + z = 0.0; + double ppos[3]; + int img[3]; + +#ifdef MASS + + double mass = 0.0; + double M; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + M = part[i].p.mass; + mass += M; + + memcpy(ppos, part[i].r.p, 3*sizeof(double)); + memcpy(img, part[i].l.i, 3*sizeof(int)); + unfold_position(ppos, img); + + x += M*ppos[0]; + y += M*ppos[1]; + z += M*ppos[2]; + } + } + + sdata[0] = x; + sdata[1] = y; + sdata[2] = z; + sdata[3] = mass; + +#else + + int npart = 0; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + npart++; + + memcpy(ppos, part[i].r.p, 3*sizeof(double)); + memcpy(img, part[i].l.i, 3*sizeof(int)); + unfold_position(ppos, img); + + x += ppos[0]; + y += ppos[1]; + z += ppos[2]; + } + } + + sdata[0] = x; + sdata[1] = y; + sdata[2] = z; + sdata[3] = (double)npart; + +#endif +} + +/* Calculate the CMS velocity of the system */ +void local_system_CMS_velocity( double *sdata ) { + int c, np, i; + Particle *part; + Cell *cell; + double x = 0.0, + y = 0.0, + z = 0.0; + +#ifdef MASS + + double mass = 0.0; + double M; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + M = part[i].p.mass; + mass += M; + + x += M*part[i].m.v[0]; + y += M*part[i].m.v[1]; + z += M*part[i].m.v[2]; + } + } + + sdata[0] = x; + sdata[1] = y; + sdata[2] = z; + sdata[3] = mass; + +#else + + int npart = 0; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + npart++; + + x += part[i].m.v[0]; + y += part[i].m.v[1]; + z += part[i].m.v[2]; + } + } + + sdata[0] = x; + sdata[1] = y; + sdata[2] = z; + sdata[3] = (double)npart; + +#endif +} + +/* Remove the CMS velocity */ +void local_galilei_transform( double *sdata ) { + int c, np, i; + Particle *part; + Cell *cell; + + for (c = 0; c < local_cells.n; c++) { + cell = local_cells.cell[c]; + part = cell->part; + np = cell->n; + + for(i = 0; i < np; i++) { + part[i].m.v[0] -= sdata[0]; + part[i].m.v[1] -= sdata[1]; + part[i].m.v[2] -= sdata[2]; + } + } +} + +#endif diff --git a/src/galilei.h b/src/galilei.h new file mode 100644 index 00000000000..c813753a67f --- /dev/null +++ b/src/galilei.h @@ -0,0 +1,52 @@ +/* + Copyright (C) 2010,2011,2012 The ESPResSo project + Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + Max-Planck-Institute for Polymer Research, Theory Group + + This file is part of ESPResSo. + + ESPResSo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ESPResSo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#ifndef GALILEI_H +#define GALILEI_H +/** \file galilei.h + * + */ + +#include "utils.h" +#include "particle_data.h" + +#ifdef GALILEI +/** broadcasts reaction parameters and sets up an entry in the ia_params, so + that the verlet radius is equal or bigger than the reaction range. +**/ +void local_kill_particle_motion(); +void local_kill_particle_forces(); +void local_system_CMS( double * ); +void local_system_CMS_velocity( double * ); +void local_galilei_transform( double * ); + +typedef struct { + + double cms[3]; + double cms_vel[3]; + +} galilei_struct; + +/** Galilei parameters. */ +extern galilei_struct gal; + +#endif + +#endif diff --git a/src/tcl/galilei_tcl.c b/src/tcl/galilei_tcl.c new file mode 100644 index 00000000000..038159b7143 --- /dev/null +++ b/src/tcl/galilei_tcl.c @@ -0,0 +1,216 @@ +/* + Copyright (C) 2010,2011,2012 The ESPResSo project + Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + Max-Planck-Institute for Polymer Research, Theory Group + + This file is part of ESPResSo. + + ESPResSo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ESPResSo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +/** \file reaction_tcl.c + * +*/ + +#include "galilei.h" +#include "parser.h" +#include "utils.h" +#include "grid.h" +#include "communication.h" +#include "integrate.h" + +/* ############### */ + +#ifdef GALILEI +/* Stop motion of the particles */ +int tcl_command_kill_particle_motion_print_usage(Tcl_Interp * interp){ + char buffer[256]; + sprintf(buffer, "Usage: kill_particle_motion\n"); + Tcl_AppendResult(interp, buffer, (char *)NULL); + return TCL_ERROR; +} + +/* Set the forces on the particles to zero */ +int tcl_command_kill_particle_forces_print_usage(Tcl_Interp * interp){ + char buffer[256]; + sprintf(buffer, "Usage: kill_particle_forces\n"); + Tcl_AppendResult(interp, buffer, (char *)NULL); + return TCL_ERROR; +} + +/* Calculate the CMS of the system */ +int tcl_command_system_CMS_print_usage(Tcl_Interp * interp){ + char buffer[256]; + sprintf(buffer, "Usage: system_CMS [folded]\n"); + Tcl_AppendResult(interp, buffer, (char *)NULL); + return TCL_ERROR; +} + +/* Calculate the CMS velocity of the system */ +int tcl_command_system_CMS_velocity_print_usage(Tcl_Interp * interp){ + char buffer[256]; + sprintf(buffer, "Usage: system_CMS_velocity\n"); + Tcl_AppendResult(interp, buffer, (char *)NULL); + return TCL_ERROR; +} + +/* Remove the CMS velocity of the system */ +int tcl_command_galilei_transform_print_usage(Tcl_Interp * interp){ + char buffer[256]; + sprintf(buffer, "Usage: galilei_transform\n"); + Tcl_AppendResult(interp, buffer, (char *)NULL); + return TCL_ERROR; +} +#endif + +/* ############### */ + +/* Stop motion of the particles */ +int tclcommand_kill_particle_motion(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ +#ifdef GALILEI + if (argc != 1 ) { + return tcl_command_kill_particle_motion_print_usage(interp); + } else { + if (ARG_IS_S_EXACT(0,"kill_particle_motion")) { + mpi_kill_particle_motion(); + return TCL_OK; + } else { + return tcl_command_kill_particle_motion_print_usage(interp); + } + } +#else + Tcl_AppendResult(interp, "GALILEI not compiled in!" ,(char *) NULL); + return (TCL_ERROR); +#endif +} + +/* Set the forces on the particles to zero */ +int tclcommand_kill_particle_forces(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ +#ifdef GALILEI + if (argc != 1 ) { + return tcl_command_kill_particle_forces_print_usage(interp); + } else { + if (ARG_IS_S_EXACT(0,"kill_particle_forces")) { + mpi_kill_particle_forces(); + return TCL_OK; + } else { + return tcl_command_kill_particle_forces_print_usage(interp); + } + } +#else + Tcl_AppendResult(interp, "GALILEI not compiled in!" ,(char *) NULL); + return (TCL_ERROR); +#endif +} + +/* Calculate the CMS of the system */ +int tclcommand_system_CMS(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ +#ifdef GALILEI + char buffer[256]; + double cmspos[3]; + int box[3]; + + if (argc != 1 && argc != 2 ) { + return tcl_command_system_CMS_print_usage(interp); + } else { + if (ARG_IS_S_EXACT(0,"system_CMS")) { + if ( argc == 2 ) { + if (ARG_IS_S_EXACT(1,"folded")) { + mpi_system_CMS(); + + memcpy(cmspos, gal.cms, 3*sizeof(double)); + box[0] = 0; box[1] = 0; box[2] = 0; + fold_position(cmspos, box); + + Tcl_PrintDouble(interp, cmspos[0], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, cmspos[1], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, cmspos[2], buffer); + Tcl_AppendResult(interp, buffer, (char *)NULL); + + return TCL_OK; + } else { + return tcl_command_system_CMS_print_usage(interp); + } + } else { + mpi_system_CMS(); + + Tcl_PrintDouble(interp, gal.cms[0], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, gal.cms[1], buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, gal.cms[2], buffer); + Tcl_AppendResult(interp, buffer, (char *)NULL); + + + return TCL_OK; + } + } else { + return tcl_command_system_CMS_print_usage(interp); + } + } +#else + Tcl_AppendResult(interp, "GALILEI not compiled in!" ,(char *) NULL); + return (TCL_ERROR); +#endif +} + +/* Calculate the CMS velocity of the system */ +int tclcommand_system_CMS_velocity(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ +#ifdef GALILEI + char buffer[256]; + + if (argc != 1 ) { + return tcl_command_system_CMS_velocity_print_usage(interp); + } else { + if (ARG_IS_S_EXACT(0,"system_CMS_velocity")) { + mpi_system_CMS_velocity(); + + Tcl_PrintDouble(interp, gal.cms_vel[0]/time_step, buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, gal.cms_vel[1]/time_step, buffer); + Tcl_AppendResult(interp, buffer, " ", (char *)NULL); + Tcl_PrintDouble(interp, gal.cms_vel[2]/time_step, buffer); + Tcl_AppendResult(interp, buffer, (char *)NULL); + + return TCL_OK; + } else { + return tcl_command_system_CMS_velocity_print_usage(interp); + } + } +#else + Tcl_AppendResult(interp, "GALILEI not compiled in!" ,(char *) NULL); + return (TCL_ERROR); +#endif +} + +/* Remove the CMS velocity of the system */ +int tclcommand_galilei_transform(ClientData data, Tcl_Interp * interp, int argc, char ** argv){ +#ifdef GALILEI + if (argc != 1 ) { + return tcl_command_galilei_transform_print_usage(interp); + } else { + if (ARG_IS_S_EXACT(0,"galilei_transform")) { + mpi_galilei_transform(); + return TCL_OK; + } else { + return tcl_command_galilei_transform_print_usage(interp); + } + } +#else + Tcl_AppendResult(interp, "GALILEI not compiled in!" ,(char *) NULL); + return (TCL_ERROR); +#endif +} diff --git a/src/tcl/galilei_tcl.h b/src/tcl/galilei_tcl.h new file mode 100644 index 00000000000..81b4c7fd688 --- /dev/null +++ b/src/tcl/galilei_tcl.h @@ -0,0 +1,45 @@ +/* + Copyright (C) 2010,2011,2012 The ESPResSo project + Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + Max-Planck-Institute for Polymer Research, Theory Group + + This file is part of ESPResSo. + + ESPResSo is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + ESPResSo is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ +#ifndef _GALILEI_TCL_H +#define _GALILEI_TCL_H +#include "parser.h" + +/** \name Exported Functions */ +/************************************************************/ +/*@{*/ + +/** tcl procedure to set up reactions +*/ + +int tclcommand_kill_particle_motion(ClientData data, Tcl_Interp *interp, + int argc, char **argv); +int tclcommand_kill_particle_forces(ClientData data, Tcl_Interp *interp, + int argc, char **argv); +int tclcommand_system_CMS(ClientData data, Tcl_Interp *interp, + int argc, char **argv); +int tclcommand_system_CMS_velocity(ClientData data, Tcl_Interp *interp, + int argc, char **argv); +int tclcommand_galilei_transform(ClientData data, Tcl_Interp *interp, + int argc, char **argv); + +/*@}*/ + +#endif