added Runge Kutta integration methods

docs/Smoldyn/SmoldynManual.tex

@@ -1518,7 +1518,7 @@ \section{Conformational spread reactions}
 \end{lstlisting}
 This will yield a warning in Smoldyn about there being multiple bimolecular reactions listed with the same reactants, but it is the right way to list these symmetric effects. In this example, the convention was followed that the latter reactant (and latter product) is the neighbor molecule, while the former reactant is the one that changes state.
 
-If a molecule has simultaneous conformational spread interactions with more than one other molecule, the simulated reaction rates may be too low; this effect is reduced to zero for short time steps and increases with longer time steps. Consider a potential reaction with two reaction channels and the probability of it happening by either channel individually is $p$. When the two channels are considered sequentially, the probability for the first happening should be p, while the probability for the second should $p/(1p)$, because it is the conditional probability of the second reaction happening, given that the first one did not happen. However, Smoldyn uses probability $p$ for all conformational spread reaction channels, which leads to a reaction rate that is too low. While this identical effect is addressed correctly for first order reactions and for state conversions of surface-bound molecules, it is not addressed for conformational spread reactions because it is nearly impossible for Smoldyn to figure out how many reaction channels are available for any particular conformational spread reaction.
+If a molecule has simultaneous conformational spread interactions with more than one other molecule, the simulated reaction rates may be too low; this effect is reduced to zero for short time steps and increases with longer time steps. Consider a potential reaction with two reaction channels and the probability of it happening by either channel individually is $p$. When the two channels are considered sequentially, the probability for the first happening should be p, while the probability for the second should $p/(1-p)$, because it is the conditional probability of the second reaction happening, given that the first one did not happen. However, Smoldyn uses probability $p$ for all conformational spread reaction channels, which leads to a reaction rate that is too low. While this identical effect is addressed correctly for first order reactions and for state conversions of surface-bound molecules, it is not addressed for conformational spread reactions because it is nearly impossible for Smoldyn to figure out how many reaction channels are available for any particular conformational spread reaction.
 
 Conformational spread reactions were tested with the configuration file confspread.txt. It simulates two reactions:
 \begin{lstlisting}[style=SSAC]

examples/S13_filaments/filament2d.txt

@@ -17,7 +17,7 @@ color red red
 
 time_start 0
 time_stop 50000
-time_step 0.01
+time_step 0.02
 
 frame_thickness 1
 
@@ -29,7 +29,7 @@ thickness 2
 polygon ve
 force_arrows 1 red
 kT 1
-dynamics Euler
+dynamics RK2
 standard_length 5
 standard_angle 0
 force_length 1
@@ -46,10 +46,11 @@ end_filament_type
 
 random_filament fil1 template 8  10 50 u
 random_filament fil2 template 49  20 50 u
-#random_filament fil3 template 49  30 50 u
+random_filament fil3 template 20  30 50 u
 #random_filament fil4 template 49  40 50 u
 #random_filament fil5 template 49  50 50 u
 
+
 cmd B pause
 #cmd E printFilament template fil1 xabf stdout
 

examples/S13_filaments/filament3.txt

@@ -17,7 +17,7 @@ color red white
 
 time_start 0
 time_stop 50000
-time_step 0.1
+time_step 0.3
 
 frame_thickness 0
 
@@ -29,7 +29,7 @@ force_arrows 10 red
 thickness 2
 polygon ve
 kT 0.1
-dynamics euler
+dynamics Euler
 standard_length 5
 standard_angle 0 0 0
 force_length 1
@@ -40,13 +40,20 @@ end_filament_type
 #start_filament template fil1
 #first_segment 0 0 0  5  0 0 0
 #add_segment 5  1 0 0
+#add_segment 5  0 1 0
+#add_segment 5  0 0 1
 #end_filament
 
 random_filament fil1 template 10  0 0 0 u u u
-random_filament fil2 template 49  20 50 1 u u u
-random_filament fil3 template 49  30 50 1 u u u
-random_filament fil4 template 49  40 50 1 u u u
-random_filament fil5 template 49  50 50 1 u u u
+random_filament fil2 template 20  20 50 1 u u u
+random_filament fil3 template 15  30 50 1 u u u
+random_filament fil4 template 30  40 50 1 u u u
+random_filament fil5 template 25  50 50 1 u u u
+
+start_filament_type template
+#standard_angle 0.3 0.3 0.1
+standard_angle 0.3 0 0
+end_filament_type
 
 #cmd B printFilament template fil1 xabyft stdout
 cmd B pause

source/Smoldyn/smoldyn.h

@@ -759,7 +759,9 @@ typedef struct latticesuperstruct
 enum FilamentDynamics
 {
     FDnone,
-    FDeuler
+    FDeuler,
+    FDRK2,
+    FDRK4
 };
 
 typedef struct segmentstruct
@@ -783,6 +785,11 @@ typedef struct filamentstruct
     int nseg;                           // number of segments
     segmentptr* segments;               // array of segments (nseg)
     double **nodes;                     // list of nodes (nseg+1)
+    double **nodes1;                    // nodes for RK dynamics (nseg+1)
+    double **nodes2;                    // nodes for RK4 dynamics (nseg+1)
+    double **nodesx;                    // old node locations (nseg+1)
+    double *roll1;											// roll for RK dynamics (nseg)
+    double *roll2;											// roll for RK dynamics (nseg)
     double **forces;										// list of forces on nodes (nseg+1)
     double *torques;										// list of segment torques (nseg)
     struct filamentstruct* frontend;    // what front attaches to