Develop

.gitignore

@@ -1,41 +1,80 @@
-*~
-*.a
+# Files that should be ignored by git for tracking
+# For explanation see the gitignore man page
+#
+# ignore all compiled python files
+*.pyc
+*.pyo
+# ignore all compiled C++ libraries and object files
+*.lo
 *.la
 *.lai
-*.lo
-*.loT
-*.o
 *.so
-*.pyo
-*.pyc
+*.Plo
+*.Po
+*.o
+*.a
+*.loT
+# ignore generated Bessel tables
+/CylindricalBesselTable.hpp
+/SphericalBesselTable.hpp
+
+# ignore files generated by pdflatex
+*.aux
+*.log
+*.pdf
+*.toc
+
+# ignore files generated by other tools such as pygmentize
+*.png
+
+# ignore data files generated by the examples
+*.out
+
+# ignore Makefiles.in generated by automake and Makefiles generated by configure
+Makefile.in
+/Makefile
+/doc/Makefile
+/binding/Makefile
+/samples/benchmark/Makefile
+/test/Makefile
+
+# all kinds of configuration files follow now
+# ignore files generated by configure 
+/config.h
+/config.log
+/config.status
+/stamp-h1
+# ignore scripts generated by configure (I think)
+/config.guess
+/config.sub
+/libtool
+/ltmain.sh
+# ignore files generated by autoconf
+/configure
+/depcomp
+/missing
+/py-compile
+/install-sh
+/m4/libtool.m4
+/m4/ltoptions.m4
+/m4/ltsugar.m4
+/m4/ltversion.m4
+/m4/lt~obsolete.m4
+/autom4te.cache/
+# ignore files generated by autoheader
+/config.h.in
+# ignore files generated by aclocal
+/aclocal.m4
+
+# ignore some other files/directories
+*~
 *.orig
 *.rej
 .libs
 .deps
 .*_history*
-m4
-Makefile
-Makefile.in
-depcomp
-config.h
-config.h.in
-ltmain.sh
-.deps
-config.status
-config.h.in
-libtool
-py-compile
-configure
-config.log
-config.guess
-config.sub
-stamp-h1
-autom4te.cache
-SphericalBesselTable.hpp
-CylindricalBesselTable.hpp
-missing
-aclocal.m4
-install-sh
+
+# Note sure why these are excluded (added by Morioshi)
 test/World_test
 test/array_helper_test
 test/filters_test
@@ -55,5 +94,6 @@ test/linear_algebra_test
 test/pointer_as_ref_test
 test/EGFRDSimulator_test
 test/geometry_test
-samples/benchmark/hardbody
-samples/*/*.out
+
+# ignore vi swp files open while editting
+*.swp

BDPropagator.hpp

@@ -10,6 +10,7 @@
 #include <boost/range/const_iterator.hpp>
 #include <boost/utility/enable_if.hpp>
 #include <boost/scoped_ptr.hpp>
+#include <boost/shared_ptr.hpp>
 #include "Defs.hpp"
 #include "generator.hpp"
 #include "exceptions.hpp"
@@ -97,15 +98,15 @@ class BDPropagator
         if (species.D() == 0.)
             return true;
 
-        position_type const displacement(drawR_free(species));
+        position_type const displacement(drawR_free(pp));
         position_type const new_pos(
             tx_.apply_boundary(
                 add(pp.second.position(), displacement)));
 
         particle_id_pair particle_to_update(
                 pp.first, particle_type(species.id(),
                     particle_shape_type(new_pos, species.radius()),
-                    species.D()));
+                    pp.second.structure_id(), species.D()));
         boost::scoped_ptr<particle_id_pair_and_distance_list> overlapped(
             tx_.check_overlap(particle_to_update.second.shape(),
                               particle_to_update.first));
@@ -158,9 +159,10 @@ class BDPropagator
     }
 
 private:
-    position_type drawR_free(species_type const& species)
+    position_type drawR_free(particle_id_pair const& pp)
     {
-        return tx_.get_structure(species.structure_id())->bd_displacement(std::sqrt(2.0 * species.D() * dt_), rng_);
+        species_type const species(tx_.get_species( pp.second.sid()) );
+        return tx_.get_structure(pp.second.structure_id())->bd_displacement(species.v() * dt_, std::sqrt(2.0 * species.D() * dt_), rng_);
     }
 
     bool attempt_reaction(particle_id_pair const& pp)
@@ -196,6 +198,7 @@ class BDPropagator
                             pp.first, particle_type(products[0],
                                 particle_shape_type(pp.second.position(),
                                                     s0.radius()),
+                                                    pp.second.structure_id(),
                                                     s0.D()));
                         boost::scoped_ptr<particle_id_pair_and_distance_list> overlapped(tx_.check_overlap(new_p.second.shape(), new_p.first));
                         if (overlapped && overlapped->size() > 0)
@@ -238,14 +241,17 @@ class BDPropagator
                                 throw propagation_error("no space");
                             }
 
-                            const Real rnd(rng_());
-                            length_type pair_distance(
-                                drawR_gbd(rnd, r01, dt_, D01));
-                            const position_type m(random_unit_vector() * pair_distance);
+                            boost::shared_ptr<structure_type> pp_structure( 
+                                tx_.get_structure( pp.second.structure_id()) );  
+
+                            position_type m( pp_structure->
+                                        dissociation_vector( rng_, r01, dt_, D01, s0.v() - s1.v() ) );
+
                             np0 = tx_.apply_boundary(pp.second.position()
-                                    + m * (s0.D() / D01));
+                                    - m * (s0.D() / D01));
                             np1 = tx_.apply_boundary(pp.second.position()
-                                    - m * (s1.D() / D01));
+                                    + m * (s1.D() / D01));
+
                             boost::scoped_ptr<particle_id_pair_and_distance_list> overlapped_s0(
                                 tx_.check_overlap(
                                     particle_shape_type(np0, s0.radius()),
@@ -268,8 +274,8 @@ class BDPropagator
 
                         tx_.remove_particle(pp.first);
                         const particle_id_pair
-                            npp0(tx_.new_particle(s0.id(), np0)),
-                            npp1(tx_.new_particle(s1.id(), np1));
+                            npp0(tx_.new_particle(s0.id(), pp.second.structure_id(), np0)),
+                            npp1(tx_.new_particle(s1.id(), pp.second.structure_id(), np1));
 
                         if (rrec_)
                         {
@@ -299,8 +305,8 @@ class BDPropagator
         }
 
         const species_type s0(tx_.get_species(pp0.second.sid())),
-                s1(tx_.get_species(pp1.second.sid()));
-        const length_type r01(s0.radius() + s1.radius());
+                           s1(tx_.get_species(pp1.second.sid()));
+        length_type r01(s0.radius() + s1.radius());
 
         const Real rnd(rng_());
         Real prob = 0;
@@ -309,7 +315,12 @@ class BDPropagator
                 i(boost::begin(rules)), e(boost::end(rules)); i != e; ++i)
         {
             reaction_rule_type const& r(*i);
-            const Real p(r.k() * dt_ / ((I_bd(r01, dt_, s0.D()) + I_bd(r01, dt_, s1.D())) * 4.0 * M_PI));
+
+            boost::shared_ptr<structure_type> pp_structure( 
+               tx_.get_structure( pp0.second.structure_id()) );  // we assume that both particles live on the same surface
+
+            const Real p(pp_structure->p_acceptance( r.k(), dt_, r01, subtract( pp1.second.position(), 
+                                                pp0.second.position() ), s0.D(), s1.D(), s0.v(), s1.v() ));
             BOOST_ASSERT(p >= 0.);
             prob += p;
             if (prob >= 1.)
@@ -330,6 +341,7 @@ class BDPropagator
                 switch (::size(products))
                 {
                 case 1:
+                // If the number or product of the reaction is one.
                     {
                         const species_id_type product(products[0]);
                         const species_type sp(tx_.get_species(product));
@@ -342,6 +354,7 @@ class BDPropagator
                                             pp1.second.position(),
                                             pp0.second.position()), s0.D())),
                                     (s0.D() + s1.D()))));
+
                         boost::scoped_ptr<particle_id_pair_and_distance_list> overlapped(
                             tx_.check_overlap(particle_shape_type(new_pos, sp.radius()),
                                               pp0.first, pp1.first));
@@ -362,7 +375,7 @@ class BDPropagator
 
                         remove_particle(pp0.first);
                         remove_particle(pp1.first);
-                        particle_id_pair npp(tx_.new_particle(product, new_pos));
+                        particle_id_pair npp(tx_.new_particle(product, pp0.second.structure_id(), new_pos));
                         if (rrec_)
                         {
                             (*rrec_)(
@@ -372,6 +385,7 @@ class BDPropagator
                         break;
                     }
                 case 0:
+                // If the two particles annihilate eachother. Just remove the particles.
                     remove_particle(pp0.first);
                     remove_particle(pp1.first);
                     break;