+ complex<double> V=0;
+
+ for (int nu=0; nu<4; nu++)
+ V += conf[ (*nb)[x*4+nu] ].phi;
+
+ const double r = exp( - M * norm(conf[x].phi)
+ - 1./M * norm(V)
+ + 2 * real( conf[x].phi * conj(V) ) );
+
+ conf[x].phi = sqrt(std::log( 1./(1-r) )) / sqrt(M)
+ * complex<double>( cos(theta), sin(theta) )
+ + V / M;