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}
 
diff --git a/myconfig-sample.h.in b/myconfig-sample.h.in
index dcd21a5b2fe..d0e42caffe1 100644
--- a/myconfig-sample.h.in
+++ b/myconfig-sample.h.in
@@ -62,6 +62,7 @@
 #define BOND_CONSTRAINT
 #define MODES
 #define BOND_VIRTUAL
+#define LANGEVIN_PER_PARTICLE
 
 /* To use address, also activate VIRTUAL_SITES and VIRTUAL_SITES_COM */
 #define ADRESS
diff --git a/src/communication.c b/src/communication.c
index 75520e2a8cd..c38a4207e05 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,75 @@ void mpi_bcast_max_mu() {
 #endif
 }
 
+#ifdef LANGEVIN_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();
+}
+#endif
+
+void mpi_set_particle_temperature_slave(int pnode, int part)
+{
+#ifdef LANGEVIN_PER_PARTICLE
+  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
+}
+
+#ifdef LANGEVIN_PER_PARTICLE
+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();
+}
+#endif
+
+void mpi_set_particle_gamma_slave(int pnode, int part)
+{
+#ifdef LANGEVIN_PER_PARTICLE
+  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();
+#endif
+}
+
 /*********************** MAIN LOOP for slaves ****************/
 
 void mpi_loop()
diff --git a/src/communication.h b/src/communication.h
index c77014df712..22538ff60d2 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 LANGEVIN_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..4a75d9f1ff8 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 LANGEVIN_PER_PARTICLE
+  part->T = -1.0;
+  part->gamma = -1.0;
+#endif
 }
 
 void free_particle(Particle *part) {
@@ -1700,6 +1705,54 @@ int tclcommand_part_parse_unfix(Tcl_Interp *interp, int argc, char **argv,
 
 #endif
 
+#ifdef LANGEVIN_PER_PARTICLE
+int part_parse_temp(Tcl_Interp *interp, int argc, char **argv,
+			 int part_num, int * change)
+{
+  double T;
+
+  *change = 1;
+
+  if (argc < 1) {
+    Tcl_AppendResult(interp, "temp 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 +2109,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 LANGEVIN_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 {
       Tcl_AppendResult(interp, "unknown particle parameter \"",
 		       argv[0],"\"", (char *)NULL);
@@ -2458,6 +2518,45 @@ int set_particle_torque(int part, double torque[3])
 
 #endif
 
+#ifdef LANGEVIN_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..ffdf763ef22 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 LANGEVIN_PER_PARTICLE
+  double T;
+  double gamma;
+#endif
+
 } Particle;
 
 /** List of particles. The particle array is resized using a sophisticated
@@ -514,6 +520,22 @@ int set_particle_dipm(int part, double dipm);
 int set_particle_virtual(int part,int isVirtual);
 #endif
 
+#ifdef LANGEVIN_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 gamma its new frictional coefficient.
+    @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.
     @param part  the particle.
diff --git a/src/thermostat.h b/src/thermostat.h
index 16f413d74b6..67d92538ca4 100644
--- a/src/thermostat.h
+++ b/src/thermostat.h
@@ -128,7 +128,10 @@ int tclcallback_thermo_ro(Tcl_Interp *interp, void *_data);
 MDINLINE void friction_thermo_langevin(Particle *p)
 {
   extern double langevin_pref1, langevin_pref2;
-
+#ifdef LANGEVIN_PER_PARTICLE
+  double langevin_pref1_temp, langevin_pref2_temp;
+#endif
+  
   int j;
 #ifdef MASS
   double massf = sqrt(PMASS(*p));
@@ -162,8 +165,29 @@ MDINLINE void friction_thermo_langevin(Particle *p)
 //    if (!(p->l.ext_flag & COORD_FIXED(j)))
     if (1==1)
 #endif
-      {
+    {
+#ifdef LANGEVIN_PER_PARTICLE  
+      if(p->gamma >= 0.) {
+        langevin_pref1_temp = -p->gamma/time_step;
+        
+        if(p->T >= 0.)
+          langevin_pref2_temp = sqrt(24.0*p->T*p->gamma/time_step);
+        else
+          langevin_pref2_temp = sqrt(24.0*temperature*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 {
+        if(p->T >= 0.)
+          langevin_pref2_temp = sqrt(24.0*p->T*langevin_gamma/time_step);
+        else          
+          langevin_pref2_temp = langevin_pref2;
+        
+        p->f.f[j] = langevin_pref1*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
-}