From 2925d6436b466a2dfed8bbe4e1aa370d5bed675f Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Fri, 25 Nov 2011 18:42:34 +0100 Subject: [PATCH 1/6] Started development of a test script for the constraint type rhomboid. Simulation works. Missing: check of results and caller interface. --- testsuite/constraint_rhomboid.tcl | 46 +++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 testsuite/constraint_rhomboid.tcl diff --git a/testsuite/constraint_rhomboid.tcl b/testsuite/constraint_rhomboid.tcl new file mode 100644 index 00000000000..3289d6964f2 --- /dev/null +++ b/testsuite/constraint_rhomboid.tcl @@ -0,0 +1,46 @@ +setmd skin 0.1 +setmd box_l 15 15 15 +setmd time_step 0.01 +setmd periodic 0 0 0 +thermostat off + +constraint rhomboid corner 5 5 5 a 5 0 0 b 0 5 0 c 0 0 5 direction outside type 0 + +part 0 pos 7.5 7.5 0 v 0 0 1 type 0 +part 1 pos 0 7.5 7.5 v 1 0 0 type 0 +part 2 pos 7.5 0 7.5 v 0 1 0 type 0 +part 3 pos 7.5 7.5 15 v 0 0 -1 type 0 +part 4 pos 15 7.5 7.5 v -1 0 0 type 0 +part 5 pos 7.5 15 7.5 v 0 -1 0 type 0 + +part 6 pos 0 0 0 v 1 1 1 type 0 +part 7 pos 0 15 0 v 1 -1 1 type 0 +part 8 pos 0 15 15 v 1 -1 -1 type 0 +part 9 pos 0 0 15 v 1 1 -1 type 0 +part 10 pos 15 0 0 v -1 1 1 type 0 +part 11 pos 15 15 0 v -1 -1 1 type 0 +part 12 pos 15 15 15 v -1 -1 -1 type 0 +part 13 pos 15 0 15 v -1 1 -1 type 0 + +part 14 pos 7.5 0 0 v 0 1 1 type 0 +part 15 pos 0 7.5 0 v 1 0 1 type 0 +part 16 pos 0 0 7.5 v 1 1 0 type 0 +part 17 pos 7.5 15 15 v 0 -1 -1 type 0 +part 18 pos 15 7.5 15 v -1 0 -1 type 0 +part 19 pos 15 15 7.5 v -1 -1 0 type 0 +part 20 pos 7.5 0 15 v 0 1 -1 type 0 +part 21 pos 0 7.5 15 v 1 0 -1 type 0 +part 22 pos 15 0 7.5 v -1 1 0 type 0 +part 23 pos 0 15 7.5 v 1 -1 0 type 0 +part 24 pos 7.5 15 0 v 0 -1 1 type 0 +part 25 pos 15 7.5 0 v -1 0 1 type 0 + +inter 0 0 lennard-jones 1. 1. 1.1225 0.25 0 + +prepare_vmd_connection "vmd" 2000 1 1 + +while {1} { + integrate 1 + imd positions + after 5 +} From a5aabe5ca8aed64537cc6b8df31a364678a25590 Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Mon, 28 Nov 2011 19:09:19 +0100 Subject: [PATCH 2/6] Stefan took and old hack of mine that allowed to scale the temperature per particle and generalized it so that one can set an arbitrary temperature and frictional coefficient for the langevin integrator for every particle. I then took this code, cleaned it up, made it switchable by a compiler switch TEMPERATURE_PER_PARTICLE and inserted it into a copy of the current master. Also made sure all demonstration scripts work with this version. To come: performance tests to decide whether to make this feature switchable or permanently on; documentation. --- src/communication.c | 66 +++++++++++++++++++++- src/communication.h | 8 +++ src/particle_data.c | 94 +++++++++++++++++++++++++++++++ src/particle_data.h | 21 +++++-- src/thermostat.h | 14 ++++- testsuite/constraint_rhomboid.tcl | 46 --------------- 6 files changed, 196 insertions(+), 53 deletions(-) delete mode 100644 testsuite/constraint_rhomboid.tcl diff --git a/src/communication.c b/src/communication.c index 75520e2a8cd..5755f282989 100644 --- a/src/communication.c +++ b/src/communication.c @@ -134,6 +134,8 @@ typedef void (SlaveCallback)(int node, int param); CB(mpi_send_vs_relative_slave) \ CB(mpi_recv_fluid_populations_slave) \ CB(mpi_recv_fluid_boundary_flag_slave) \ + CB(mpi_set_particle_temperature_slave) \ + CB(mpi_set_particle_gamma_slave) \ // create the forward declarations #define CB(name) void name(int node, int param); @@ -2518,7 +2520,6 @@ void mpi_iccp3m_iteration_slave(int dummy, int dummy2) #endif } - /********************* REQ_ICCP3M_INIT********/ int mpi_iccp3m_init(int n_induced_charges) { @@ -2600,6 +2601,69 @@ void mpi_bcast_max_mu() { #endif } +#ifdef TEMPERATURE_PER_PARTICLE +/******************** REQ_SEND_PARTICLE_T ********************/ +void mpi_set_particle_temperature(int pnode, int part, double _T) +{ + mpi_call(mpi_set_particle_temperature_slave, pnode, part); //TODO: really? + + if (pnode == this_node) { + Particle *p = local_particles[part]; + /* here the setting actually happens, if the particle belongs to the local node */ + p->T = _T; + } + else { + MPI_Send(&_T, 1, MPI_DOUBLE, pnode, SOME_TAG, MPI_COMM_WORLD); + } + + on_particle_change(); +} + +void mpi_set_particle_temperature_slave(int pnode, int part) +{ + double s_buf = 0.; + if (pnode == this_node) { + Particle *p = local_particles[part]; + MPI_Status status; + MPI_Recv(&s_buf, 1, MPI_DOUBLE, 0, SOME_TAG, MPI_COMM_WORLD, &status); + /* here the setting happens for nonlocal nodes */ + p->T = s_buf; + } + + on_particle_change(); +} +#endif + +void mpi_set_particle_gamma(int pnode, int part, double gamma) +{ + mpi_call(mpi_set_particle_gamma_slave, pnode, part); + + if (pnode == this_node) { + Particle *p = local_particles[part]; + /* here the setting actually happens, if the particle belongs to the local node */ + p->gamma = gamma; + } + else { + MPI_Send(&gamma, 1, MPI_DOUBLE, pnode, SOME_TAG, MPI_COMM_WORLD); + } + + on_particle_change(); +} + +void mpi_set_particle_gamma_slave(int pnode, int part) +{ + double s_buf = 0.; + if (pnode == this_node) { + Particle *p = local_particles[part]; + MPI_Status status; + MPI_Recv(&s_buf, 1, MPI_DOUBLE, 0, SOME_TAG, MPI_COMM_WORLD, &status); + /* here the setting happens for nonlocal nodes */ + p->gamma = s_buf; + } + + on_particle_change(); +} + /*********************** MAIN LOOP for slaves ****************/ void mpi_loop() diff --git a/src/communication.h b/src/communication.h index c77014df712..cd9ff07c664 100644 --- a/src/communication.h +++ b/src/communication.h @@ -384,6 +384,14 @@ void mpi_bcast_coulomb_params(); /** Issue REQ_SEND_EXT: send nex external flag and external force. */ void mpi_send_ext(int pnode, int part, int flag, int mask, double force[3]); +#ifdef TEMPERATURE_PER_PARTICLE +/** Issue REQ_SEND_PARTICLE_T: send particle type specific temperature. */ +void mpi_set_particle_temperature(int pnode, int part, double _T); + +/** Issue REQ_SEND_PARTICLE_T: send particle type specific frictional coefficient. */ +void mpi_set_particle_gamma(int pnode, int part, double gamma); +#endif + /** Issue REQ_BCAST_COULOMB: send new coulomb parameters. */ void mpi_bcast_constraint(int del_num); diff --git a/src/particle_data.c b/src/particle_data.c index cfbec123791..f439c67121e 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -1700,6 +1700,54 @@ int tclcommand_part_parse_unfix(Tcl_Interp *interp, int argc, char **argv, #endif +#ifdef TEMPERATURE_PER_PARTICLE +int part_parse_T(Tcl_Interp *interp, int argc, char **argv, + int part_num, int * change) +{ + double T; + + *change = 1; + + if (argc < 1) { + Tcl_AppendResult(interp, "T_scaling requires 1 argument", (char *) NULL); + return TCL_ERROR; + } + /* set temperature */ + if (! ARG_IS_D(0, T)) + return TCL_ERROR; + + if (set_particle_temperature(part_num, T) == TCL_ERROR) { + Tcl_AppendResult(interp, "set particle position first", (char *)NULL); + return TCL_ERROR; + } + + return TCL_OK; +} + +int part_parse_gamma(Tcl_Interp *interp, int argc, char **argv, + int part_num, int * change) +{ + double gamma; + + *change = 1; + + if (argc < 1) { + Tcl_AppendResult(interp, "gamma requires 1 argument", (char *) NULL); + return TCL_ERROR; + } + /* set temperature scaling factor */ + if (! ARG_IS_D(0, gamma)) + return TCL_ERROR; + + if (set_particle_gamma(part_num, gamma) == TCL_ERROR) { + Tcl_AppendResult(interp, "set particle position first", (char *)NULL); + return TCL_ERROR; + } + + return TCL_OK; +} +#endif + int tclcommand_part_parse_bond(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) { @@ -2056,6 +2104,13 @@ int tclcommand_part_parse_cmd(Tcl_Interp *interp, int argc, char **argv, err = tclcommand_part_parse_exclusion(interp, argc-1, argv+1, part_num, &change); #endif +#ifdef TEMPERATURE_PER_PARTICLE + else if (ARG0_IS_S("temp")) + err = part_parse_T(interp, argc-1, argv+1, part_num, &change); + else if (ARG0_IS_S("gamma")) + err = part_parse_gamma(interp, argc-1, argv+1, part_num, &change); +#endif + else { Tcl_AppendResult(interp, "unknown particle parameter \"", argv[0],"\"", (char *)NULL); @@ -2458,6 +2513,45 @@ int set_particle_torque(int part, double torque[3]) #endif +#ifdef TEMPERATURE_PER_PARTICLE +int set_particle_temperature(int part, double T) +{ + 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_set_particle_temperature(pnode, part, T); + return TCL_OK; +} + +int set_particle_gamma(int part, double gamma) +{ + 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_set_particle_gamma(pnode, part, gamma); + return TCL_OK; +} +#endif + #ifdef EXTERNAL_FORCES int set_particle_ext(int part, int flag, double force[3]) { diff --git a/src/particle_data.h b/src/particle_data.h index 0cf20cad462..84598847a49 100644 --- a/src/particle_data.h +++ b/src/particle_data.h @@ -229,6 +229,12 @@ typedef struct { /** list of particles, with which this particle has no nonbonded interactions */ IntList el; #endif + +#ifdef TEMPERATURE_PER_PARTICLE + double T; + double gamma; +#endif + } Particle; /** List of particles. The particle array is resized using a sophisticated @@ -505,13 +511,20 @@ int set_particle_dip(int part, double dip[3]); int set_particle_dipm(int part, double dipm); #endif -#ifdef VIRTUAL_SITES -/** Call only on the master node: set particle dipole moment (absolut value). +#ifdef TEMPERATURE_PER_PARTICLE +/** Call only on the master node: set particle temperature. + @param part the particle. + @param T its new temperature. + @return TCL_OK if particle existed +*/ +int set_particle_temperature(int part, double T); + +/** Call only on the master node: set particle frictional coefficient. @param part the particle. - @param isVirtual its new dipole moment. + @param gamma its new frictional coefficient. @return TCL_OK if particle existed */ -int set_particle_virtual(int part,int isVirtual); +int set_particle_gamma(int part, double gamma); #endif #ifdef EXTERNAL_FORCES diff --git a/src/thermostat.h b/src/thermostat.h index 16f413d74b6..edecaf807a1 100644 --- a/src/thermostat.h +++ b/src/thermostat.h @@ -127,8 +127,12 @@ int tclcallback_thermo_ro(Tcl_Interp *interp, void *_data); */ MDINLINE void friction_thermo_langevin(Particle *p) { +#ifdef TEMPERATURE_PER_PARTICLE + double langevin_pref1_temp, langevin_pref2_temp; //TODO: Maybe unneeded. +#else extern double langevin_pref1, langevin_pref2; - +#endif + int j; #ifdef MASS double massf = sqrt(PMASS(*p)); @@ -162,8 +166,14 @@ MDINLINE void friction_thermo_langevin(Particle *p) // if (!(p->l.ext_flag & COORD_FIXED(j))) if (1==1) #endif - { + { +#ifdef TEMPERATURE_PER_PARTICLE + langevin_pref1_temp = -p->gamma/time_step; + langevin_pref2_temp = sqrt(24.0*p->T*p->gamma/time_step); + p->f.f[j] = langevin_pref1_temp*p->m.v[j]*PMASS(*p) + langevin_pref2_temp*(d_random()-0.5)*massf; +#else p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf; +#endif } #ifdef EXTERNAL_FORCES else p->f.f[j] = 0; diff --git a/testsuite/constraint_rhomboid.tcl b/testsuite/constraint_rhomboid.tcl deleted file mode 100644 index 3289d6964f2..00000000000 --- a/testsuite/constraint_rhomboid.tcl +++ /dev/null @@ -1,46 +0,0 @@ -setmd skin 0.1 -setmd box_l 15 15 15 -setmd time_step 0.01 -setmd periodic 0 0 0 -thermostat off - -constraint rhomboid corner 5 5 5 a 5 0 0 b 0 5 0 c 0 0 5 direction outside type 0 - -part 0 pos 7.5 7.5 0 v 0 0 1 type 0 -part 1 pos 0 7.5 7.5 v 1 0 0 type 0 -part 2 pos 7.5 0 7.5 v 0 1 0 type 0 -part 3 pos 7.5 7.5 15 v 0 0 -1 type 0 -part 4 pos 15 7.5 7.5 v -1 0 0 type 0 -part 5 pos 7.5 15 7.5 v 0 -1 0 type 0 - -part 6 pos 0 0 0 v 1 1 1 type 0 -part 7 pos 0 15 0 v 1 -1 1 type 0 -part 8 pos 0 15 15 v 1 -1 -1 type 0 -part 9 pos 0 0 15 v 1 1 -1 type 0 -part 10 pos 15 0 0 v -1 1 1 type 0 -part 11 pos 15 15 0 v -1 -1 1 type 0 -part 12 pos 15 15 15 v -1 -1 -1 type 0 -part 13 pos 15 0 15 v -1 1 -1 type 0 - -part 14 pos 7.5 0 0 v 0 1 1 type 0 -part 15 pos 0 7.5 0 v 1 0 1 type 0 -part 16 pos 0 0 7.5 v 1 1 0 type 0 -part 17 pos 7.5 15 15 v 0 -1 -1 type 0 -part 18 pos 15 7.5 15 v -1 0 -1 type 0 -part 19 pos 15 15 7.5 v -1 -1 0 type 0 -part 20 pos 7.5 0 15 v 0 1 -1 type 0 -part 21 pos 0 7.5 15 v 1 0 -1 type 0 -part 22 pos 15 0 7.5 v -1 1 0 type 0 -part 23 pos 0 15 7.5 v 1 -1 0 type 0 -part 24 pos 7.5 15 0 v 0 -1 1 type 0 -part 25 pos 15 7.5 0 v -1 0 1 type 0 - -inter 0 0 lennard-jones 1. 1. 1.1225 0.25 0 - -prepare_vmd_connection "vmd" 2000 1 1 - -while {1} { - integrate 1 - imd positions - after 5 -} From cdf72a04f3ea4a83c2a3739dd3915619ac755dc7 Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Mon, 28 Nov 2011 19:23:19 +0100 Subject: [PATCH 3/6] Readded accidentaly overwritten code segment for virtual sites in 'particle_data.h' --- src/particle_data.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/particle_data.h b/src/particle_data.h index 84598847a49..53df3887bc2 100644 --- a/src/particle_data.h +++ b/src/particle_data.h @@ -511,6 +511,15 @@ int set_particle_dip(int part, double dip[3]); int set_particle_dipm(int part, double dipm); #endif +#ifdef VIRTUAL_SITES +/** Call only on the master node: set particle dipole moment (absolut value). + @param part the particle. + @param isVirtual its new dipole moment. + @return TCL_OK if particle existed +*/ +int set_particle_virtual(int part,int isVirtual); +#endif + #ifdef TEMPERATURE_PER_PARTICLE /** Call only on the master node: set particle temperature. @param part the particle. From 4097351426dc512d9db640876f06e9d78ee24abf Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Tue, 29 Nov 2011 17:15:27 +0100 Subject: [PATCH 4/6] Initialize particles with invalid T and gamma. Langevin integrator decides, based on this information, whether to take the global thermostat parameters or the particle specific ones for thermalization. --- src/communication.c | 8 +++++++- src/particle_data.c | 17 +++++++++++------ src/thermostat.h | 17 ++++++++++------- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/communication.c b/src/communication.c index 5755f282989..122d8a3397f 100644 --- a/src/communication.c +++ b/src/communication.c @@ -2618,9 +2618,11 @@ void mpi_set_particle_temperature(int pnode, int part, double _T) on_particle_change(); } +#endif void mpi_set_particle_temperature_slave(int pnode, int part) { +#ifdef TEMPERATURE_PER_PARTICLE double s_buf = 0.; if (pnode == this_node) { Particle *p = local_particles[part]; @@ -2631,9 +2633,10 @@ void mpi_set_particle_temperature_slave(int pnode, int part) } on_particle_change(); -} #endif +} +#ifdef TEMPERATURE_PER_PARTICLE void mpi_set_particle_gamma(int pnode, int part, double gamma) { mpi_call(mpi_set_particle_gamma_slave, pnode, part); @@ -2649,9 +2652,11 @@ void mpi_set_particle_gamma(int pnode, int part, double gamma) on_particle_change(); } +#endif void mpi_set_particle_gamma_slave(int pnode, int part) { +#ifdef TEMPERATURE_PER_PARTICLE double s_buf = 0.; if (pnode == this_node) { Particle *p = local_particles[part]; @@ -2662,6 +2667,7 @@ void mpi_set_particle_gamma_slave(int pnode, int part) } on_particle_change(); +#endif } /*********************** MAIN LOOP for slaves ****************/ diff --git a/src/particle_data.c b/src/particle_data.c index f439c67121e..fd675a0f5f8 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -94,7 +94,7 @@ void auto_exclusion(int distance); * particle initialization functions ************************************************/ -void init_particle(Particle *part) +void init_particle(Particle *part) { /* ParticleProperties */ part->p.identity = -1; @@ -206,6 +206,11 @@ void init_particle(Particle *part) #ifdef ADRESS part->p.adress_weight = 1.0; #endif + +#ifdef TEMPERATURE_PER_PARTICLE + part->T = -1.0; + part->gamma = -1.0; +#endif } void free_particle(Particle *part) { @@ -1701,7 +1706,7 @@ int tclcommand_part_parse_unfix(Tcl_Interp *interp, int argc, char **argv, #endif #ifdef TEMPERATURE_PER_PARTICLE -int part_parse_T(Tcl_Interp *interp, int argc, char **argv, +int part_parse_temp(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) { double T; @@ -1709,7 +1714,7 @@ int part_parse_T(Tcl_Interp *interp, int argc, char **argv, *change = 1; if (argc < 1) { - Tcl_AppendResult(interp, "T_scaling requires 1 argument", (char *) NULL); + Tcl_AppendResult(interp, "temp requires 1 argument", (char *) NULL); return TCL_ERROR; } /* set temperature */ @@ -2106,9 +2111,9 @@ int tclcommand_part_parse_cmd(Tcl_Interp *interp, int argc, char **argv, #ifdef TEMPERATURE_PER_PARTICLE else if (ARG0_IS_S("temp")) - err = part_parse_T(interp, argc-1, argv+1, part_num, &change); - else if (ARG0_IS_S("gamma")) - err = part_parse_gamma(interp, argc-1, argv+1, part_num, &change); + err = part_parse_temp(interp, argc-1, argv+1, part_num, &change); + else if (ARG0_IS_S("gamma")) + err = part_parse_gamma(interp, argc-1, argv+1, part_num, &change); #endif else { diff --git a/src/thermostat.h b/src/thermostat.h index edecaf807a1..e768d8f1619 100644 --- a/src/thermostat.h +++ b/src/thermostat.h @@ -127,10 +127,9 @@ int tclcallback_thermo_ro(Tcl_Interp *interp, void *_data); */ MDINLINE void friction_thermo_langevin(Particle *p) { -#ifdef TEMPERATURE_PER_PARTICLE - double langevin_pref1_temp, langevin_pref2_temp; //TODO: Maybe unneeded. -#else extern double langevin_pref1, langevin_pref2; +#ifdef TEMPERATURE_PER_PARTICLE + double langevin_pref1_temp, langevin_pref2_temp; #endif int j; @@ -167,10 +166,14 @@ MDINLINE void friction_thermo_langevin(Particle *p) if (1==1) #endif { -#ifdef TEMPERATURE_PER_PARTICLE - langevin_pref1_temp = -p->gamma/time_step; - langevin_pref2_temp = sqrt(24.0*p->T*p->gamma/time_step); - p->f.f[j] = langevin_pref1_temp*p->m.v[j]*PMASS(*p) + langevin_pref2_temp*(d_random()-0.5)*massf; +#ifdef TEMPERATURE_PER_PARTICLE + if(p->gamma >= 0. && p->T >= 0.) { + langevin_pref1_temp = -p->gamma/time_step; + langevin_pref2_temp = sqrt(24.0*p->T*p->gamma/time_step); + p->f.f[j] = langevin_pref1_temp*p->m.v[j]*PMASS(*p) + langevin_pref2_temp*(d_random()-0.5)*massf; + } + else + p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf; #else p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf; #endif From eb3f52d034cf3353ab9849f1200abb21fc5fda8f Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Tue, 29 Nov 2011 21:28:30 +0100 Subject: [PATCH 5/6] Deleted compiler switch 'TEMPERATURE_PER_PARTICLE' after thorough investigation of the performance impact of the feature. The result was that for 10-1000 ideal gas particles, the slowdown is (0.8+/-1.7)%. --- src/communication.c | 8 -------- src/communication.h | 2 -- src/particle_data.c | 13 ++----------- src/particle_data.h | 5 ----- src/thermostat.h | 6 ------ 5 files changed, 2 insertions(+), 32 deletions(-) diff --git a/src/communication.c b/src/communication.c index 122d8a3397f..08900c48fbe 100644 --- a/src/communication.c +++ b/src/communication.c @@ -2601,7 +2601,6 @@ void mpi_bcast_max_mu() { #endif } -#ifdef TEMPERATURE_PER_PARTICLE /******************** REQ_SEND_PARTICLE_T ********************/ void mpi_set_particle_temperature(int pnode, int part, double _T) { @@ -2618,11 +2617,9 @@ void mpi_set_particle_temperature(int pnode, int part, double _T) on_particle_change(); } -#endif void mpi_set_particle_temperature_slave(int pnode, int part) { -#ifdef TEMPERATURE_PER_PARTICLE double s_buf = 0.; if (pnode == this_node) { Particle *p = local_particles[part]; @@ -2633,10 +2630,8 @@ void mpi_set_particle_temperature_slave(int pnode, int part) } on_particle_change(); -#endif } -#ifdef TEMPERATURE_PER_PARTICLE void mpi_set_particle_gamma(int pnode, int part, double gamma) { mpi_call(mpi_set_particle_gamma_slave, pnode, part); @@ -2652,11 +2647,9 @@ void mpi_set_particle_gamma(int pnode, int part, double gamma) on_particle_change(); } -#endif void mpi_set_particle_gamma_slave(int pnode, int part) { -#ifdef TEMPERATURE_PER_PARTICLE double s_buf = 0.; if (pnode == this_node) { Particle *p = local_particles[part]; @@ -2667,7 +2660,6 @@ void mpi_set_particle_gamma_slave(int pnode, int part) } on_particle_change(); -#endif } /*********************** MAIN LOOP for slaves ****************/ diff --git a/src/communication.h b/src/communication.h index cd9ff07c664..e0855c501d2 100644 --- a/src/communication.h +++ b/src/communication.h @@ -384,13 +384,11 @@ void mpi_bcast_coulomb_params(); /** Issue REQ_SEND_EXT: send nex external flag and external force. */ void mpi_send_ext(int pnode, int part, int flag, int mask, double force[3]); -#ifdef TEMPERATURE_PER_PARTICLE /** Issue REQ_SEND_PARTICLE_T: send particle type specific temperature. */ void mpi_set_particle_temperature(int pnode, int part, double _T); /** Issue REQ_SEND_PARTICLE_T: send particle type specific frictional coefficient. */ void mpi_set_particle_gamma(int pnode, int part, double gamma); -#endif /** Issue REQ_BCAST_COULOMB: send new coulomb parameters. */ void mpi_bcast_constraint(int del_num); diff --git a/src/particle_data.c b/src/particle_data.c index fd675a0f5f8..1f364a4b503 100644 --- a/src/particle_data.c +++ b/src/particle_data.c @@ -207,10 +207,8 @@ void init_particle(Particle *part) part->p.adress_weight = 1.0; #endif -#ifdef TEMPERATURE_PER_PARTICLE part->T = -1.0; part->gamma = -1.0; -#endif } void free_particle(Particle *part) { @@ -1705,7 +1703,6 @@ int tclcommand_part_parse_unfix(Tcl_Interp *interp, int argc, char **argv, #endif -#ifdef TEMPERATURE_PER_PARTICLE int part_parse_temp(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) { @@ -1751,7 +1748,6 @@ int part_parse_gamma(Tcl_Interp *interp, int argc, char **argv, return TCL_OK; } -#endif int tclcommand_part_parse_bond(Tcl_Interp *interp, int argc, char **argv, int part_num, int * change) @@ -2109,13 +2105,10 @@ int tclcommand_part_parse_cmd(Tcl_Interp *interp, int argc, char **argv, err = tclcommand_part_parse_exclusion(interp, argc-1, argv+1, part_num, &change); #endif -#ifdef TEMPERATURE_PER_PARTICLE else if (ARG0_IS_S("temp")) err = part_parse_temp(interp, argc-1, argv+1, part_num, &change); - else if (ARG0_IS_S("gamma")) - err = part_parse_gamma(interp, argc-1, argv+1, part_num, &change); -#endif - + else if (ARG0_IS_S("gamma")) + err = part_parse_gamma(interp, argc-1, argv+1, part_num, &change); else { Tcl_AppendResult(interp, "unknown particle parameter \"", argv[0],"\"", (char *)NULL); @@ -2518,7 +2511,6 @@ int set_particle_torque(int part, double torque[3]) #endif -#ifdef TEMPERATURE_PER_PARTICLE int set_particle_temperature(int part, double T) { int pnode; @@ -2555,7 +2547,6 @@ int set_particle_gamma(int part, double gamma) mpi_set_particle_gamma(pnode, part, gamma); return TCL_OK; } -#endif #ifdef EXTERNAL_FORCES int set_particle_ext(int part, int flag, double force[3]) diff --git a/src/particle_data.h b/src/particle_data.h index 53df3887bc2..251387ecdf6 100644 --- a/src/particle_data.h +++ b/src/particle_data.h @@ -230,11 +230,8 @@ typedef struct { IntList el; #endif -#ifdef TEMPERATURE_PER_PARTICLE double T; double gamma; -#endif - } Particle; /** List of particles. The particle array is resized using a sophisticated @@ -520,7 +517,6 @@ int set_particle_dipm(int part, double dipm); int set_particle_virtual(int part,int isVirtual); #endif -#ifdef TEMPERATURE_PER_PARTICLE /** Call only on the master node: set particle temperature. @param part the particle. @param T its new temperature. @@ -534,7 +530,6 @@ int set_particle_temperature(int part, double T); @return TCL_OK if particle existed */ int set_particle_gamma(int part, double gamma); -#endif #ifdef EXTERNAL_FORCES /** Call only on the master node: set particle external forced. diff --git a/src/thermostat.h b/src/thermostat.h index e768d8f1619..750fb408ede 100644 --- a/src/thermostat.h +++ b/src/thermostat.h @@ -128,9 +128,7 @@ int tclcallback_thermo_ro(Tcl_Interp *interp, void *_data); MDINLINE void friction_thermo_langevin(Particle *p) { extern double langevin_pref1, langevin_pref2; -#ifdef TEMPERATURE_PER_PARTICLE double langevin_pref1_temp, langevin_pref2_temp; -#endif int j; #ifdef MASS @@ -166,7 +164,6 @@ MDINLINE void friction_thermo_langevin(Particle *p) if (1==1) #endif { -#ifdef TEMPERATURE_PER_PARTICLE if(p->gamma >= 0. && p->T >= 0.) { langevin_pref1_temp = -p->gamma/time_step; langevin_pref2_temp = sqrt(24.0*p->T*p->gamma/time_step); @@ -174,9 +171,6 @@ MDINLINE void friction_thermo_langevin(Particle *p) } else p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf; -#else - p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) + langevin_pref2*(d_random()-0.5)*massf; -#endif } #ifdef EXTERNAL_FORCES else p->f.f[j] = 0; From 1d2cc5187c00b85cfa85ffdbfbadd837daa65eb3 Mon Sep 17 00:00:00 2001 From: Georg Rempfer Date: Wed, 30 Nov 2011 00:49:40 +0100 Subject: [PATCH 6/6] Documented the new feature 'part pid temp T gamma g' in the users guide and cross referenced the documentation on the 'part' command with the section on the Langevin thermostat. --- doc/ug/introduction.tex | 2 +- doc/ug/part.tex | 4 ++++ doc/ug/setup.tex | 4 ++++ 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/doc/ug/introduction.tex b/doc/ug/introduction.tex index 7117341e6d8..0a457e31f3a 100644 --- a/doc/ug/introduction.tex +++ b/doc/ug/introduction.tex @@ -140,7 +140,7 @@ \section{Available simulation methods} \section{Basic program structure} \label{sec:structure} -As already mentioned, \es contains of two components. +As already mentioned, \es consists of two components. The simulation engine is written in C for the sake of computational efficiency. The steering or control level is interfaced to the kernel via an interpreter diff --git a/doc/ug/part.tex b/doc/ug/part.tex index 6fab6f41883..54a541422cc 100644 --- a/doc/ug/part.tex +++ b/doc/ug/part.tex @@ -34,6 +34,8 @@ \subsection{Defining particle properties} \opt{v \var{vx} \var{vy} \var{vz}} \opt{f \var{fx} \var{fy} \var{fz}} \opt{bond \var{bondid} \var{pid2} \dots} + \opt{temp \var{T}} + \opt{gamma \var{g}} \require{1}{\opt{q \var{charge}}} \require{2}{\opt{quat \var{q1} \var{q2} \var{q3} \var{q4}}} \require{2}{\opt{omega \var{x} \var{y} \var{z}}} @@ -85,6 +87,8 @@ \subsection{Defining particle properties} be an existing particle. The \var{bondid} is used for the inter command to define bonded interactions. \item[bond delete] Will delete all bonds attached to this particle. +\item[\opt{temp \var{T}}] If used in combination with the Langevin thermostat (as documented in section \ref{sec:thermostat}), sets the temperature \var{T} individually for the particle with id \var{pid}. This allows to simulate systems containing particles of different temperatures. Caution: this has no influence on any other thermostat then the Langevin thermostat. +\item[\opt{gamma \var{g}}] If used in combination with the Langevin thermostat (as documented in section \ref{sec:thermostat}), sets the frictional coefficient \var{T} individually for the particle with id \var{pid}. This allows to simulate systems containing particles with different diffusion constants. Caution: this has no influence on any other thermostat then the Langevin thermostat. \item[\opt{q \var{charge}}] Sets the charge of this particle to $q$. \item[\opt{quat \var{q1} \var{q2} \var{q3} \var{q4}}] Sets the quaternion representation of the rotational position of this diff --git a/doc/ug/setup.tex b/doc/ug/setup.tex index 6f983f292a4..4de4ba25de6 100644 --- a/doc/ug/setup.tex +++ b/doc/ug/setup.tex @@ -180,6 +180,10 @@ \subsection{Langevin thermostat} If the feature \feature{ROTATION} is compiled in, the rotational degrees of freedom are also coupled to the thermostat. +Using the Langevin thermostat, it is posible to set a temperature and a friction +coefficient for every particle individually. Consult the \texttt{part} commands reference +(chapter \ref{chap:part}) for information on how to achieve this. + \subsection{Dissipative Particle Dynamics (DPD) } \label{sec:DPD} \index{DPD|mainindex}