diff --git a/docs/Smoldyn/SmoldynCodeDoc.pdf b/docs/Smoldyn/SmoldynCodeDoc.pdf index bce64c08..369072d4 100644 Binary files a/docs/Smoldyn/SmoldynCodeDoc.pdf and b/docs/Smoldyn/SmoldynCodeDoc.pdf differ diff --git a/docs/Smoldyn/SmoldynCodeDoc.tex b/docs/Smoldyn/SmoldynCodeDoc.tex index 66eac56e..ebeb0438 100644 --- a/docs/Smoldyn/SmoldynCodeDoc.tex +++ b/docs/Smoldyn/SmoldynCodeDoc.tex @@ -6,6 +6,7 @@ \usepackage{listings} \usepackage{url} \usepackage{hyperref} +\usepackage{amsmath} % settings for listings package \lstset{ @@ -3953,170 +3954,167 @@ \subsection*{Functions} % Filaments (functions in smolfilament.c) \section{Filaments (functions in smolfilament.c)} -Filaments are 1-dimensional structures that are designed to represent cytoskeletal filaments, such as actin or microtubules, or DNA. They are intended to be much longer than they are wide. I wrote most of the code here, although some of it was extended from code that Edward Rolls wrote; he was a graduate student in Radek Erban's lab during 2016, when I was at the Newton Institute in Cambridge. +Filaments are 1-dimensional structures that are designed to represent cytoskeletal filaments, such as actin or microtubules, or DNA. They are intended to be much longer than they are wide. Essentially all of the functions here are designed for filaments in 3D space. However, they also work for 2D filaments, provided that the filament parameters are set up so that the filament always stays within the 2D plane. \subsection*{Data structures declared in smoldyn.h} -Filaments are defined in a heirarchical fashion. At the top level, there is one filament superstructure which contains a list of filament types that are used in the simulation. Below that are multiple filament type structures. Each individual filament of a particular type looks and behaves identically. For example, actin in the cytoplasm might be one type; another might be double-stranded DNA, and yet another is single-stranded DNA. Each filament type contains multiple individual filaments. - -Each filament is composed of either segments or beads, depending on the dynamics that were chosen for that type. Both types of filaments can have stretching energies. Bead filaments never have bending energies but segment filaments can. Segment thicknesses can vary over the length of a filament, but not bead radii. Each segment or bead is only part of a single filament; it's possible for a filament to branch off of a segment or bead in another filament, but that bead is not shared but owned by the original filament. The following table lists the dynamics and beads vs. segment representations. - -\begin{longtable}[c]{ll} -Dynamics & Beads or segements\\ -\hline -rigidbeads & beads\\ -rigidsegments & segments \\ -rouse & beads\\ -alberts & segments\\ -nedelec & don't know -\end{longtable} +Filaments are defined in a heirarchical fashion. At the top level, there is one filament superstructure which contains a list of filament types that are used in the simulation. Below that are multiple filament type structures. For example, actin in the cytoplasm might be one type; another might be double-stranded DNA, and yet another is single-stranded DNA. Each filament type contains multiple individual filaments, which is the next level down, each of which looks and behaves identically. Each filament is composed of a list of segments. Segments can vary in lengths and thicknesses. \subsubsection*{Enumerations} \begin{lstlisting} -enum FilamentDynamics {FDnone, FDrouse, FDalberts, FDnedelec}; enum FilamentBiology {FBactin, FBmicrotubule, FBintermediate, FBdsDNA, FBssDNA, FBother}; +enum FilamentDynamics {FDnone, FDrouse, FDalberts, FDnedelec}; \end{lstlisting} -These enumerations are used in the filament types. The filament dynamics enumeration describes how the filament is to represented and what algorithms to use for its dynamics. The filament biology enumeration is used to describe what kind of filament this is, which can then be used to pre-populate many of the filament parameters. +The filament biology enumeration is used to describe what kind of filament this is, which can then be used to pre-populate many of the filament parameters. The filament dynamics enumeration describes what algorithms to use for its dynamics. They are summarized below; this list will be expanded. -\subsubsection*{Beads and segments} - -\begin{lstlisting} -typedef struct beadstruct { - double xyz[3]; // bead coordinates - double xyzold[3]; // bead coordinates for prior time - } *beadptr; -\end{lstlisting} +\begin{longtable}[c]{ll} +Dynamics & Description\\ +\hline +none & filament is rigid and does not move\\ +newton & use Newton ODE integration\\ +\end{longtable} -Each bead, for filaments that have them, is quite simple. A bead has its current position, in \ttt{xyz} and its old position, in \ttt{xyzold}. For 2D simulations, only the first two vales are used in each position vector. +\subsubsection*{Segments} \begin{lstlisting} typedef struct segmentstruct { - double xyzfront[3]; // Coords. for segment front - double xyzback[3]; // Coords. for segment back - double len; // segment length - double ypr[3]; // relative ypr angles - double dcm[9]; // relative dcm - double adcm[9]; // absolute segment orientation - double thk; // thickness of segment - } *segmentptr; + struct filamentstruct* fil; // owning filament + int index; // self index along filament + double xyzfront[3]; // Coords. for segment front + double xyzback[3]; // Coords. for segment back + double len; // segment length + double thk; // thickness of segment + double ypr[3]; // relative ypr angles + double dcm[9]; // relative dcm + double adcm[9]; // absolute segment orientation + } * segmentptr; \end{lstlisting} -Segments are a bit more complicated than beads. \ttt{xyzfront} is the $x, y, z$ coordinate of the segment starting point and \ttt{xyzback} is the $x, y, z$ coordinate of the segment end point. \ttt{len} is the segment length. \ttt{ypr} is the relative yaw, pitch, and roll angle for the segment relative to the orientation of the prior segment. \ttt{dcm} is the relative direction cosine matrix for the segment relative to the orientation of the one before it. Identical information is contained in \ttt{ypr} and \ttt{dcm}, but the latter is here for faster computation. \ttt{adcm} is the direction cosine matrix for the absolute orientation. \ttt{thk} is the thickness of the segment. For 2D simulations, only the first 2 values are used in the coordinates vectors and only the first value is used in the ypr angles. +Segments start by listing the filament that owns them. Then they have their own position within the filament in \ttt{index}; note that the segment array within a filament is shifted by \ttt{frontseg}, so a pointer to the same segment would be found with \ttt{segment->fil->segments[segment->index+segment->fil->frontseg]} and the result would be identical to \ttt{segment}. + +\ttt{xyzfront} is the $x, y, z$ coordinate of the segment starting point and \ttt{xyzback} is the $x, y, z$ coordinate of the segment end point. \ttt{len} is the segment length and \ttt{thk} is the segment thickness. \ttt{ypr} is the relative yaw, pitch, and roll angle for the segment relative to the orientation of the prior segment. \ttt{dcm} is the relative direction cosine matrix for the segment relative to the orientation of the one before it. Identical information is contained in \ttt{ypr} and \ttt{dcm}, but the latter is here for faster computation. \ttt{adcm} is the direction cosine matrix for the absolute orientation; again, this information is redundant but included for faster computation. For 2D simulations, only the first 2 values are used in the coordinates vectors and the third coordinate should always equal 0; only the first value (yaw) is used in the ypr angles and the other two should always equal 0; and only the first 2 rows and columns are used in the direction cosine matrices while the others equal 0 for unequal indices and 1 for equal indices. -At present, filaments are designed for all segments of the same length. This wasn't what I had in mind initially, but may be sufficiently easier to program that it's worth keeping this constraint.\\ +Segments lengths can vary but standard lengths are constant over an entire filament. This means that segments all have roughly equal lengths. I initially thought that it could be useful to have large standard length variations to focus computation on more important regions, but this is sufficiently hard to program that it's not worth it.\\ \subsubsection*{Filaments} \begin{lstlisting} typedef struct filamentstruct { - struct filamenttypestruct *filtype; // filament type structure - char *filname; // filament name (reference, not owned) - int maxbs; // number of beads or segments allocated - int nbs; // number of beads or segments - int frontbs; // index of front bead or segment - int backbs; // index of back bead or segment - beadptr *beads; // array of beads if any - segmentptr *segments; // array of segments if any - struct filamentstruct *frontend; // what front attaches to if anything - struct filamentstruct *backend; // what back attaches to if anything - int maxbranch; // max number of branches off this filament - int nbranch; // number of branches off this filament - int *branchspots; // list of bead or segments where branches are - struct filamentstruct **branches; // list of branching filaments - int maxmonomer; // number of monomers allocated - int nmonomer; // number of monomers - int frontmonomer; // index of front monomer - int backmonomer; // index of back monomer - char *monomers; // monomer codes - } *filamentptr; + struct filamenttypestruct* filtype; // owning filament type + char* filname; // filament name (ref, not owned) + int maxseg; // number of segments allocated + int nseg; // number of segments + int frontseg; // index of front segment + segmentptr* segments; // array of segments if any + struct filamentstruct* frontend; // what front attaches to + struct filamentstruct* backend; // what back attaches to + int maxbranch; // max branches off this filament + int nbranch; // num branches off this filament + int* branchspots; // segments where branches are + struct filamentstruct** branches; // list of branching filaments + int maxmonomer; // number of monomers allocated + int nmonomer; // number of monomers + int frontmonomer; // index of front monomer + char* monomers; // monomer codes + } * filamentptr; \end{lstlisting} -Each filament data structure represents a single unbranched chain of either beads or segments. Multiple filaments of the same type (e.g. two microtubules) are stored in multiple data structures. In the filament data structure, \ttt{filtype} points to the filament type that owns this filament, \ttt{filindex} tells what the index number is of this particular filament, and \ttt{filname} points to the name of this filament; names are assigned by default if not assigned by the user. +Each filament data structure represents a single unbranched chain of segments. Multiple filaments of the same type (e.g. two microtubules) are stored in multiple data structures. In the filament data structure, \ttt{filtype} points to the filament type that owns this filament and \ttt{filname} points to the name of this filament; names are assigned by default if not assigned by the user. -Whether a filament is composed of beads or segments depends on its dynamics type, and is given in the \ttt{isbead} element of the type. Either way, \ttt{maxbs} and \ttt{nbs} give the allocated and actual numbers of either beads or segments. These are stored in an array, where the index of the front is in \ttt{frontbs} and the index of the back is in \ttt{backbs-1}. These are implemented so that the back index is always greater than the front index; I considered using a circular array, but this is easier to use and runs faster in most algorithms. If the filament is composed of beads, then they are stored in \ttt{beads} and if it is composed of segments, they are stored in \ttt{segments}. +Next, \ttt{maxseg} and \ttt{nseg} give the allocated and actual numbers of segments. These are stored in an array, where the index of the front is in \ttt{frontseg}. This is not a circular array, which is more complicated to use, but allows for the front index to be non-zero because this leads to simpler and faster operation if there is regular assembly and/or disassembly. Segments are stored in \ttt{segments}. -If \ttt{backbs==frontbs}, then this is an empty filament with no beads or segments at all. If \ttt{backbs==frontbs+1}, then this is a single element filament; it has one bead if it is bead style, or one segment if it is segment style. +Filaments are allowed to branch. If the current filament is a branch off another filament, then \ttt{frontend} points to the filament that the front end attaches to and/or \ttt{backend} points to the filament that the back end attaches to; these are \ttt{NULL} if they are free ends. Additionally, this filament knows about the filaments that branch off this one, which are listed in \ttt{maxbranch}, \ttt{nbranch}, and \ttt{branches}. The branched filaments branch off at the locations listed in \ttt{branchspots}. -Filaments are allowed to branch. If the current filament is a branch off of another filament, then \ttt{frontend} points to the filament that the front end attaches to and/or \ttt{backend} points to the filament that the back end attaches to; these are \ttt{NULL} if they are free ends. Additionally, this filament knows about the filaments that branch off of this one, which are listed in \ttt{maxbranch}, \ttt{nbranch}, and \ttt{branches}. The branched filaments branch off at the locations listed in \ttt{branchspots}. - -Along the length of each filament are equally spaced monomer codes, which represent things like DNA sequence or modification state of proteins. They are stored in a non-circular array with parameters \ttt{maxmonomer}, \ttt{nmonomer}, \ttt{frontmonomer}, \ttt{backmonomer}, and then the list \ttt{monomers}. This array is analogous to the beads or segments array, but may be larger or smaller because there's no necessary correlation between monomer codes and beads or segments.\\ +Along the length of each filament are equally spaced monomer codes, which represent things like DNA sequence or modification state of proteins. They are stored in an array with parameters \ttt{maxmonomer}, \ttt{nmonomer}, \ttt{frontmonomer}, \ttt{backmonomer}, and then the list \ttt{monomers}. This array is analogous to the segments array, but may be larger or smaller because there's no necessary correlation between monomer codes and segments.\\ \subsubsection*{Filament types} \begin{lstlisting} typedef struct filamenttypestruct { - struct filamentsuperstruct *filss; // filament superstructure - char *ftname; // filament type name - enum FilamentDynamics dynamics; // dynamics for the filament - int isbead; // 1 for bead model, 0 for segments - enum FilamentBiology biology; // Biological name for filament type - int bundlevalue; // number of microfilaments in bundle - double color[4]; // filament color - double edgepts; // thickness of edge for drawing - unsigned int edgestipple[2]; // edge stippling [factor, pattern] - enum DrawMode drawmode; // polygon drawing mode - double shiny; // shininess - double stdlen; // minimum energy segment length - double stdypr[3]; // minimum energy bend angle - double klen; // force constant for length - double kypr[3]; // force constant for angle - int maxface; // number of filament faces allocated - int nface; // number of filament faces - char **facename; // list of face names - double facetwist; // twisting rate of faces along filament - double kT; // thermodynamic temperature, [0, inf) - double treadrate; // treadmilling rate constant - double viscosity; // viscosity - double beadradius; // bead radius - int maxfil; // maximum number of filaments - int nfil; // actual number of filaments - filamentptr *fillist; // list of filaments - char **filnames; // names of filaments -} *filamenttypeptr; + struct filamentsuperstruct* filss; // filament superstructure + char* ftname; // fil type name (ref, not owned) + enum FilamentDynamics dynamics; // dynamics for the filament + enum FilamentBiology biology; // biological name for filaments + double bundlevalue; // microfilaments in bundle + double color[4]; // filament color + double edgepts; // thickness of edge for drawing + unsigned int edgestipple[2]; // edge stippling [factor, pattern] + enum DrawMode drawmode; // polygon drawing mode + double shiny; // shininess + double stdlen; // minimum energy segment length + double stdypr[3]; // minimum energy bend angle + double klen; // force constant for length + double kypr[3]; // force constant for angle + double kT; // thermodynamic temp, [0,inf) + double treadrate; // treadmilling rate constant + double viscosity; // viscosity + double filradius; // segment radius + int maxface; // filament faces allocated + int nface; // number of filament faces + char** facename; // list of face names + double facetwist; // twisting rate of faces + int maxfil; // maximum number of filaments + int nfil; // actual number of filaments + int autonamenum; // next number for autoname + filamentptr* fillist; // list of filaments + char** filnames; // names of filaments + } * filamenttypeptr; \end{lstlisting} Filament types describe various types of filaments, such as actin, microtubule, etc. All filaments of the same type have the same graphical display, force constants, dynamic simulation methods, etc. It's reasonable to have multiple types that are represent a single biological type. For example, there could be actin in cytoplasm, actin in nucleoplasm, and an actin bundle. -The \ttt{filss} element points to the filament superstructure that owns this filament type, and the name of the type is in \ttt{ftname}. The type of simulation dynamics is in \ttt{dynamics}, where the options are listed above. Whether these dynamics work for beads or segments is fixed by the type, but, for convenience, it is copied over into \ttt{isbead} so that routines can easily determine if this is a bead or segment representation. The \ttt{biology} element is not needed, but is useful because a user can specifiy the biology and then Smoldyn auto-fills the rest of the parameters. The \ttt{bundlevalue} parameter gives the number of filaments that are bundled together in this particular type. +The \ttt{filss} element points to the filament superstructure that owns this filament type, and the name of the type is in \ttt{ftname}. The type of simulation dynamics is in \ttt{dynamics}, where the options are listed above. The \ttt{biology} element is not needed, but is useful because a user can specifiy the biology and then Smoldyn auto-fills the rest of the parameters. The \ttt{bundlevalue} parameter gives the number of microfilaments that are bundled together in this particular type. In the graphical display parameters, \ttt{color} is the filament color, with red, green, blue, and alpha values. \ttt{edgepts} is the edge thickness for drawing, \ttt{edgestipple} is the stippling code for the edge stippleing, if any. \ttt{drawmode} is the polymer drawing mode, which can be face, edge, vertex, of a combination of these. \ttt{shiny} is the shininess for graphical display. These graphical display statements are very similar to those used for surfaces. -Filament mechanics information is in the next elements. \ttt{stdlen} is the standard segment length, meaning the segment length at minimum energy, where segments are neither stretched nor compressed. \ttt{stdypr} is the minimum energy relative ypr angle. \ttt{klen} is the stretching force constant. Its value is $< 0$ for an infinite force constant, meaning that the segment lengths are fixed. \ttt{kypr} is the bending force constant. For any element, its value is 0 for zero force constant, meaning a freely jointed chain, and $< 0$ for an infinite force constant, meaning a fixed bending angle. \ttt{kT} is the thermodynamic energy. This parameter isn't ideal here because it allows a different temperature for different components in the simulation, but it's here for now. \ttt{treadrate} is the treadmilling rate constant. \ttt{viscosity} is the medium viscosity for computing drag coefficients. \ttt{filradius} is the radius of either beads or segments, depending on the representation. For 2D filaments, \ttt{stdypr} values 1 and 2 should have value 0 and \ttt{kypr} values 1 and 2 should have value -1 (i.e. bending out of the $x, y$ plane can't happen, and all that's left is the yaw option). +Filament mechanics information is in the next elements. \ttt{stdlen} is the standard segment length, meaning the segment length at minimum energy, where segments are neither stretched nor compressed. \ttt{stdypr} is the minimum energy relative ypr angle. \ttt{klen} is the stretching force constant. Its value is $< 0$ for an infinite force constant, meaning that the segment lengths are fixed. \ttt{kypr} is the bending force constant. For any element, its value is 0 for zero force constant, meaning a freely jointed chain, and $< 0$ for an infinite force constant, meaning a fixed bending angle. \ttt{kT} is the thermodynamic energy. This parameter isn't ideal here because it allows a different temperature for different components in the simulation, but it's here for now. \ttt{treadrate} is the treadmilling rate constant. \ttt{viscosity} is the medium viscosity for computing drag coefficients. \ttt{filradius} is the radius of segments. -Finally, \ttt{maxfil} and \ttt{nfil} give the number of allocated and actual filaments of this type, and they are listed in \ttt{fillist}. The filaments have names in \ttt{filnames}.\\ +Next, \ttt{maxfil} and \ttt{nfil} give the number of allocated and actual filaments of this type, and they are listed in \ttt{fillist}. The filaments have names in \ttt{filnames}. Names can be assigned automatically, in which case they are named as the filament type name appended with an integer; this integer is found from \ttt{autonamenum}.\\ \subsubsection*{Filament superstructure} \begin{lstlisting} typedef struct filamentsuperstruct { - enum StructCond condition; // structure condition - struct simstruct *sim; // simulation structure - int maxtype; // maximum number of filament types - int ntype; // actual number of filament types - char **ftnames; // filament type names - filamenttypeptr *filtypes; // list of filament types - } *filamentssptr; + enum StructCond condition; // structure condition + struct simstruct* sim; // simulation structure + int maxtype; // maximum number of filament types + int ntype; // actual number of filament types + char** ftnames; // filament type names + filamenttypeptr* filtypes; // list of filament types + } * filamentssptr; \end{lstlisting} The superstructure contains a list of filament types. It starts with its condition element in \ttt{condition} and a pointer up the heirarchy to the simulation structure in \ttt{sim}. The list of filament types is allocated with \ttt{maxtype} entries, has \ttt{ntype} entries, and has list \ttt{ftnames} for the type names and \ttt{filtypes} is the actual list of types.\\ \subsection*{Filament math} -Following are the basic equations for the filament relative angles ($\mathbf{A}$, \ttt{dcm}), absolute angles ($\mathbf{B}$, \ttt{adcm}), and positions ($\mathbf{x}$, \ttt{xyz}). See Andrews, 2014 for a much more complete description of filament math. +Following are the basic equations for the filament relative angles ($\mathbf{A}$, \ttt{dcm}), absolute angles ($\mathbf{B}$, \ttt{adcm}), and positions ($\mathbf{x}$, \ttt{xyz}). See also Goldstein, Andrews, 2004 and Andrews, 2014. + +Consider direction cosine matrix $\mathbf{\Phi}$, which can be expanded as +\begin{equation} +\mathbf{\Phi} = \left[ \begin{array}{ccc} \Phi_{Xx} & \Phi_{Yx} & \Phi_{Zx} \\ \Phi_{Xy} & \Phi_{Yy} & \Phi_{Zy} \\ \Phi_{Xz} & \Phi_{Yz} & \Phi_{Zz} \end{array} \right] +\label{eq:DirectionCosineMatrix} +\end{equation} +One interpretation of the direction cosine matrix is that it converts between two frames of reference. Here, $X$, $Y$, and $Z$ are the unit vectors of the ``lab frame'', meaning the absolute coordinates of the simulation, and $x$, $y$, and $z$ are the unit vectors of the ``molecule frame'', meaning the local coordinates of a filament segment. The direction cosine matrix expresses the dot products of these unit vectors. As a result, each row of the matrix gives the lab frame coordinates of each unit vector of the molecule frame, and each column of the matrix gives the molecule frame coordinates of each unit vector of the lab frame. It is unitary, meaning that all its eigenvalues equal 1, and its transpose is its inverse, so $\mathbf{\Phi}^T \mathbf{\Phi} = \mathbf{1}$, where $\mathbf{1}$ is the identity matrix. Suppose $\mathbf{B}$ is a direction cosine matrix (called an absolute one here) and $\mathbf{p}$ is a column vector in the molecule frame. For example, $\mathbf{p}^T=[1,0,0]$ represents the direction parallel to the segment. Also, $\mathbf{p}_L$ is the same vector but expressed in the lab frame. These vectors can be converted back and forth with +\begin{align} +\mathbf{p}_L^T &= \mathbf{p}^T \cdot \mathbf{B} \hspace{2cm} \mathbf{p}^T = \mathbf{p}_L^T \cdot \mathbf{B}^T \\ +\mathbf{p}_L &= \mathbf{B}^T \cdot \mathbf{p} \hspace{2cm} \mathbf{p} = \mathbf{B}\cdot \mathbf{p}_L +\end{align} +Note that the Sphere.c library has the function \ttt{Sph\_DcmtxUnit}. -\begin{math} -\mathbf{B}_0 = \mathbf{A}_0 -\end{math} +A different interpretation of the direction cosine matrix is that it represents a rotation of a vector in a constant coordinate system, in which it can be called a rotation matrix. Refering back to eq. \ref{eq:DirectionCosineMatrix}, $X$, $Y$, and $Z$ represent the vector direction before rotation and $x$, $y$, and $z$ represent the vector after rotation. Left-multiplying the rotation matrix on a column vector yields the rotated version of the vector. -\begin{math} -\mathbf{B}_i = \mathbf{A}_i \cdot \mathbf{B}_{i-1} -\end{math} +For a filament, an unrotated version of segment 0, meaning that it's aligned with the lab coordinate system, lies along direction $[1,0,0]^T$. This gets rotated into its actual direction using the absolute direction cosine matrix $\mathbf{B}_0$, and $\mathbf{A}_0$ is defined to be the same, +$$\mathbf{B}_0 = \mathbf{A}_0$$ +Right-multiplying either of these, which are obviously the same, by $[1,0,0]^T$ gives the actual segment direction in the lab coordinate system. Segment 1 is rotated from segment 0 by the rotation matrix $\mathbf{A}_1$. Its absolute direction cosine matrix is +$$\mathbf{B}_1 = \mathbf{A}_1 \cdot \mathbf{B}_0$$ +Again, right-multiplying by $[1,0,0]^T$ gives the segment direction in the lab coordinate system. This pattern continues with additional segments, so +$$\mathbf{B}_i = \mathbf{A}_i \cdot \mathbf{B}_{i-1} = \mathbf{A}_i \mathbf{A}_{i-1} \cdots \mathbf{A}_1 \mathbf{B}_0$$ +We can also rotate the other direction, going from the back of the filament toward the front, with +$$\mathbf{B}_i = \mathbf{A}^T_{i+1} \cdot \mathbf{B}_{i+1}.$$ -\begin{math} -\mathbf{B}_i = \mathbf{A}^T_{i+1} \cdot \mathbf{B}_{i+1} -\end{math} \subsection*{Function declarations.} As much as possible, functions are declared locally rather than in the smoldynfuncs.h header file. This simplifies the code reading because it clarifies which functions might be called externally versus those that are only called internally. @@ -4125,6 +4123,12 @@ \subsection*{Function declarations.} \item[\underline{Enumerations}] +\item[\ttt{char *filFB2string(enum FilamentBiology fb, char *string)}] +Local. Converts filament biology type to string. Returns ``none'' for unrecognized biology type. Writes result to \ttt{string} and also returns it directly. + +\item[\ttt{enum FilamentBiology filstring2FB(const char *string)}] +Local. Converts filament biology string to enumerated biology type. Returns \ttt{FBnone} for unrecognized input. + \item[\ttt{char *filFD2string(enum FilamentDynamics fd, char *string);}] \hfill \\ Local. Converts filament dynamics type to string. Returns ``none'' for unrecognized dynamics type. Writes result to \ttt{string} and also returns it directly. @@ -4133,57 +4137,44 @@ \subsection*{Function declarations.} \hfill \\ Local. Converts filament dynamics string to enumerated dynamics type. Returns \ttt{FDnone} for unrecognized input. -\item[\ttt{char *filFB2string(enum FilamentBiology fb, char *string)}] -Local. Converts filament biology type to string. Returns ``other'' for unrecognized biology type. Writes result to \ttt{string} and also returns it directly. - -\item[\ttt{enum FilamentBiology filstring2FB(const char *string)}] -Local. Converts filament biology string to enumerated biology type. Returns \ttt{FBother} for unrecognized input. - \item[\underline{Low level utilities}] \item[\ttt{double}] \ttt{filRandomLength(const filamenttypeptr filtype, double thickness, double sigmamult);} \hfill \\ -Local. Returns a random segment length using the mechanics parameters given in filament type \ttt{filtype}. \ttt{thickness} is the thickness of the new segment and \ttt{sigmamult} is multiplied by the normal standard deviation of the length. If $k_{len} \leq 0$, this returns the standard length; otherwise, the returned length is Gaussian distributed, with mean equal to \ttt{fil->lstd} and standard deviation, which is constrained to be always positive, equal to +Local. Returns a random segment length using the mechanics parameters given in filament type \ttt{filtype}. \ttt{thickness} is the thickness of the new segment and \ttt{sigmamult} is multiplied by the normal standard deviation of the length. If $k_{len} < 0$, this returns the standard length; otherwise, the returned length is Gaussian distributed, with mean equal to \ttt{fil->lstd} and standard deviation, which is constrained to be always positive, equal to $$\sigma_{mult}\sqrt{\frac{kT}{R*k_{len}}}$$ \item[\ttt{double}] -\ttt{*filRandomAngle(const filamenttypeptr filtype, double thickness, double *angle, double sigmamult);} +\ttt{*filRandomAngle(const filamenttypeptr filtype,int dim,int n,double thickness,double sigmamult,double *angle));} \hfill \\ -Local. Returns a random bending angle, which is a relative ypr angle, in \ttt{angle} and directly for filament type \ttt{filtype}, without changing the filament. The new segment has thickness \ttt{thickness} and the computed standard deviation is multiplied by \ttt{sigmamult}. If $k_{ypr}$ is less than 0 for a particular coordinate, this returns the standard bending angle; if it is equal to 0, this returns a uniformly distributed bending angle; and if it is greater than 0, this returns a mean of $ypr_{std}$ and standard deviation of +Local. Returns a random bending angle, which is a relative ypr angle, in \ttt{angle} (which needs to be allocated with 3 elements) and directly. This is for filament type \ttt{filtype}, system dimensionality \ttt{dim} and segment number \ttt{n}; the number is only checked to see if it's zero, in which case angles are chosen uniformly over all possibilities or non-zero, in which case angles are chosen from angle bending forces. Angles are for a segment with thickness \ttt{thickness} and the computed standard deviation is multiplied by \ttt{sigmamult}. For \ttt{n>0}, if $k_{ypr}$ is less than 0 for a particular coordinate, this returns the standard bending angle, if it is equal to 0, this returns a uniformly distributed bending angle, and if it is greater than 0, this returns a mean of $ypr_{std}$ and standard deviation of $$\sigma_{mult} \sqrt{\frac{kT}{R*k_{ypr}}}$$ +For 2D space, this sets the last two elements of \ttt{angle} to 0. \item[\underline{Computations on filaments}] \item[\ttt{double filStretchEnergy(const filamentptr fil, int seg1, int seg2);}] \hfill \\ -Computes the stretching energy for filament \ttt{fil}. Enter \ttt{seg1} as the first segment to compute from, or as -1 to indicate that it should start at the front of the filament, and \ttt{set2} as the last segment to compute to, or as -1 to indicate that calculation should end at the end of the filament. The equation is (the thickness only applies to segment models): +Computes the stretching energy for filament \ttt{fil}. Enter \ttt{seg1} as the first segment to compute from, or as -1 to indicate that it should start at the front of the filament, and \ttt{seg2} as one more than the last segment to compute to, or as -1 to indicate that calculation should end at the end of the filament. The equation is $$E_{stretch} = \sum_{s} \frac{R_s k_{len} (l_s - l_{std})^2}{2}$$ \item[\ttt{double filBendEnergy(const filamentptr fil, int seg1, int seg2);}] \hfill \\ -Computes the bending energy for filament \ttt{fil}. This returns 0 for bead filaments because they always have 0 bending energy. Enter \ttt{seg1} as the first segment to compute from, or as -1 to indicate that it should start at the front of the filament, and \ttt{seg2} as the last segment to compute to, or as -1 to indicate that calculation should end at the end of the filament. The equation is +Computes the bending energy for filament \ttt{fil}. Enter \ttt{seg1} as the first segment to compute from, or as -1 to indicate that it should start at the front of the filament, and \ttt{seg2} as the last segment to compute to, or as -1 to indicate that calculation should end at the end of the filament. Here, \ttt{seg2} is inclusive of the last segment; the bends that are included here start at the back of \ttt{seg1} and continue up to the front of \ttt{seg2}; there are one fewer bends in a filament than the number of segments. The equation is $$E_{bend}=\sum_{s=1}^{n} \frac{R_{s-1}+R_s}{2} \frac{k_y(a_y-a_{y, std})^2+k_p(a_p-a_{p, std})^2+k_r(a_r-a_{r, std})^2}{2}$$ -The first term in the sum computes the average thickness of the segment in front of and behind the bend. The other term is the squared bending angle on each coordintate. +The first term in the sum computes the average thickness of the segment in front of and behind the bend. The other term is the squared bending angle on each coordinate (for 3D; only the first term is included for 2D). \item[\underline{Memory management}] -\item[\ttt{beadptr beadalloc();}] -\hfill \\ -Allocates memory for a single bead and initializes this bead. Returns a pointer to this bead, or \ttt{NULL} if memory could not be allocated. - -\item[\ttt{void beadfree(beadptr bead);}] -\hfill \\ -Frees memory for a single bead. - \item[\ttt{segmentptr segmentalloc();}] \hfill \\ Allocates memory for a single segment and initializes this segment. Returns a pointer to this segment, or \ttt{NULL} if memory could not be allocated. The segment length is initialized to 0 to indicate that it hasn't been set up yet, and the thickness is initialized to 1, which is good default value. Ypr angles are initialized to 0 and direction cosine matrices to unit matrices. @@ -4193,9 +4184,9 @@ \subsection*{Function declarations.} Frees memory for a single segment. \item[\ttt{filamentptr}] -\ttt{filalloc(filamentptr fil, int isbead, int maxbs, int maxbranch, int maxmonomer);} +\ttt{filalloc(filamentptr fil, int maxseg, int maxbranch, int maxmonomer);} \hfill \\ -Allocates or expands a filament, returning a pointer to the result. This does not shrink any elements of filaments. Enter \ttt{isbead} as 1 if this is a bead type filament and 0 if it is segment type (this only matters for allocating beads or segments), \ttt{maxbs} to expand the number of beads or segments, \ttt{maxbranch} to expand the number of branches off of this filament, and \ttt{maxmonomer} to expand the monomer list. This returns \ttt{NULL} if memory could not be allocated but does not free things, which can cause a minor memory leak, but is better than risking a segfault and easier than a sophisticated analysis of what can and cannot be freed. +Allocates or expands a filament, returning a pointer to the result. This does not shrink any elements of filaments. Enter \ttt{maxseg} to expand the number of segments, \ttt{maxbranch} to expand the number of branches off of this filament, and \ttt{maxmonomer} to expand the monomer list. This returns \ttt{NULL} if memory could not be allocated but does not free things, which can cause a minor memory leak, but is better than risking a segfault and easier than a sophisticated analysis of what can and cannot be freed. \item[\ttt{void filfree(filamentptr fil);}] \hfill \\ @@ -4204,7 +4195,7 @@ \subsection*{Function declarations.} \item[\ttt{filamenttypeptr}] \ttt{filamenttypealloc(filamenttypeptr filtype, int maxfil, int maxface);} \hfill \\ -Allocates or expands memory for a filament type structure and intilizes it. Returns a pointer to this filament type, or \ttt{NULL} if memory could not be allocated. Enter \ttt{filtype} as an existing one if it is to be expanded, \ttt{maxfil} as the number of filament spaces to be allocated, and \ttt{maxface} as the number of face names to be allocated. This only expands the filament or face lists if the requested sizes are larger than the current sizes. When putting new filaments into the list, they are allocated here and initialized, but have no segments or beads. +Allocates or expands memory for a filament type structure and intilizes it. Returns a pointer to this filament type, or \ttt{NULL} if memory could not be allocated. Enter \ttt{filtype} as an existing one if it is to be expanded, \ttt{maxfil} as the number of filament spaces to be allocated, and \ttt{maxface} as the number of face names to be allocated. This only expands the filament or face lists if the requested sizes are larger than the current sizes. When putting new filaments into the list, they are allocated here and initialized, but have no segments. \item[\ttt{void filamenttypefree(filamenttypeptr filtype);}] \hfill \\ @@ -4223,10 +4214,10 @@ \subsection*{Function declarations.} \item[\ttt{void filoutput(const filamentptr fil);}] \hfill \\ -Outputs all of the key information about a filament to the display. +Outputs all of the key information about a filament to the display, including essential information about each segment. \item[\ttt{void filtypeoutput(const filamenttypeptr filtype);}] -Outputs all of the key information about a filament type to the display. +Outputs all of the key information about a filament type to the display, along with information about each filament. \item[\ttt{void filssoutput(const simptr sim);}] \hfill \\ @@ -4234,27 +4225,30 @@ \subsection*{Function declarations.} \item[\ttt{void filwrite(const simptr sim, FILE *fptr);}] \hfill \\ -Writes filament information to a file in Smoldyn format for loading in again later on. This function isn't written yet. +Writes filament information to a file in Smoldyn format for loading in again later on. NOT WRITTEN YET. \item[\ttt{int filcheckparams(const simptr sim, int *warnptr);}] \hfill \\ -Checks filament parameters to make sure they are all reasonable. This function isn't written yet. +Checks filament parameters to make sure they are all reasonable. NOT WRITTEN YET. \item[\underline{Filament manipulation}] -\item[\ttt{void filArrayShift(filamentptr fil, int shift);}] +\item[\ttt{void filUpdateSegmentIndices(filamentptr fil);}] \hfill \\ -Shifts the bead or segment list in the filament (beads if they are defined and segments otherwise) either to higher indices or to lower indices, adjusting the \ttt{frontbs} and \ttt{backbs} elements to account for the shift. Enter \ttt{shift} as a positive number to increase all of the indices, as a negative number to decrease all of the indices, and as 0 to have the filament centered in the allocated memory. If \ttt{shift} is made too large or too small, such that the filament would move off the end of the allocated memory, then it is automatically fixed so that the filament just gets shifted to the end of the allocated memory instead. Thus, an easy way to make sure that the filament is at the end of the allocated memory is to enter \ttt{shift} as $\pm$\ttt{fil->maxbs}. +Overwrites the segment \ttt{index} element for each segment in a filament to make them correct. Call this, if needed, after adding or removing segments to or from a filament. -\item[\ttt{int}] -\ttt{filAddBead(filamentptr fil, const double *x, double length, char endchar)} +\item[\ttt{int filGetFilIndex(simptr sim,const char *name,int *ftptr);}] \hfill \\ -Adds a bead to filament \ttt{fil}. If this is the first bead, then \ttt{x} needs to be set to the starting location of the filament; otherwise, \ttt{x} is ignored. \ttt{length} is the length of the distance between beads, and \ttt{endchar} should be set to `b' to add to the back of the filament or `f' to add to the front of the filament. Relative angles are randomly chosen from a circle for 2D or a sphere for 3D, to create a freely-jointed chain. +Returns the index of the filament named \ttt{name} or -1 if not found; also, returns the filament type in \ttt{ftptr}. This searches through all filament types for the filament with this name. If multiple filaments with the same name are found, this returns -2. \ttt{ftptr} should be sent in pointing to a memory address, but the contents are unimportant. Example: \ttt{filGetFilIndex(sim,"fil1",\&ft);}. + +\item[\ttt{void filArrayShift(filamentptr fil, int shift);}] +\hfill \\ +Shifts the segment list in the filament either to higher indices or to lower indices, adjusting the \ttt{frontseg} element to account for the shift. Enter \ttt{shift} as a positive number to increase all the indices, as a negative number to decrease all the indices, and as 0 to have the filament centered in the allocated memory. If \ttt{shift} is made too large or too small, such that the filament would move off the end of the allocated memory, then it is automatically fixed so that the filament just gets shifted to the end of the allocated memory instead. Thus, an easy way to make sure that the filament is at the end of the allocated memory is to enter \ttt{shift} as $\pm$\ttt{fil->maxseg}. \item[\ttt{int}] \ttt{filAddSegment(filamentptr fil, const double *x, double length, const double *angle, double thickness, char endchar);} \hfill \\ -Adds a segment to filament \ttt{fil}. If this is the first segment, then \ttt{x} needs to be set to the starting location of the filament; otherwise, \ttt{x} is ignored. \ttt{length} is the length of the segment, \ttt{angle} is the relative angle of the segment facing along the filament from front to back (for the first segment, this is the absolute angle), \ttt{thickness} is the segment thickness, and \ttt{endchar} should be set to `b' to add to the back of the filament or `f' to add to the front of the filament. For segments added to the front, \ttt{angle} is the new angle from the new first monomer to the next monomer, facing towards the back. If the first segment is added to the back, then \ttt{angle} is the angle of the new segment off of the coordinate system. If the first segment is added to the front, then \ttt{angle} is the angle from the new segment to the coordinate system. +Adds a segment to filament \ttt{fil}. If this is the first segment, then \ttt{x} becomes the starting location of the filament; otherwise, \ttt{x} is ignored. \ttt{length} is the length of the segment, \ttt{angle} is the relative angle of the segment facing along the filament from front to back (for the first segment, this is the absolute angle), \ttt{thickness} is the segment thickness, and \ttt{endchar} should be set to `b' to add to the back of the filament or `f' to add to the front of the filament. For segments added to the front, \ttt{angle} is the new angle from the new first monomer to the next monomer, facing towards the back. If the first segment is added to the back, then \ttt{angle} is the angle of the new segment off of the coordinate system. If the first segment is added to the front, then \ttt{angle} is the angle from the new segment to the coordinate system. The \ttt{angle} vector must have 3 elements in it, even for 2D simulations; in the latter case, the latter two elements should be set to 0. Returns 0 for success and 1 for out of memory. For the first segment: $\mathbf{x}_0=\ttt{x}$, @@ -4275,48 +4269,65 @@ \subsection*{Function declarations.} \item[\ttt{int}] \ttt{filAddRandomSegments(filamentptr fil, int number, const char *xstr, const char *ystr, const char *zstr, const char *thtstr, const char* phistr, const char* chistr, double thickness)} \hfill \\ -Add \ttt{number} of random segments to filament \ttt{fil}. Enter \ttt{xstr}, \ttt{ystr}, and \ttt{zstr} as strings for the starting position coordinates of a new filament, and \ttt{thtstr}, \ttt{phistr}, and \ttt{chistr} as string for the starting angle of a new filament. Here, math equations are allowed, using Smoldyn variables, or use the `u' character to indicate uniform starting position within the system volume or angle in spherical coordinates. \ttt{thickness} is the thickness of the segments being added. +Add \ttt{number} of random segments to the back filament \ttt{fil}, while ignoring any surfaces or other potential collisions. Enter \ttt{xstr}, \ttt{ystr}, and \ttt{zstr} as strings for the starting position coordinates of a new filament, and \ttt{thtstr}, \ttt{phistr}, and \ttt{chistr} as strings for the starting angle of a new filament. Here, math equations are allowed, using Smoldyn variables, or use the `u' character to indicate uniform starting position within the system volume or angle in spherical coordinates. \ttt{thickness} is the thickness of the segments being added. For 2D simulations, the $z$, $\phi$, and $\chi$ values are automatically set to zero, regardless of what is entered in \ttt{zstr}, \ttt{phistr}, and \ttt{chistr}. -Returns 0 for success, 2 for error in positions, or 3 for error in angles. - -\item[\ttt{int}] -\ttt{filAddRandomBeads(filamentptr fil, int number, const char *xstr, const char *ystr, const char *zstr);} -\hfill \\ -This is the same as \ttt{filAddRandomSegments}, but adds beads instead of segments. +Returns 0 for success, 2 for syntax error in position strings, or 3 for syntax error in angle strings. \item[\ttt{int filRemoveSegment(filamentptr fil, char endchar);}] \hfill \\ -Removes one segment or one bead from either the front or back end of a filament. Specify the end in \ttt{endchar} with `f' for front and `b' for back. +Removes one segment or one bead from either the front or back end of a filament. Specify the end in \ttt{endchar} with `f' for front and `b' for back. Returns 0 for success or -1 for a filament that already had no segments. This does not shift the segments within their array. \item[\ttt{void filTranslate(filamentptr fil, const double *vect, char func)}] \hfill \\ -Translates an entire filament. Enter \ttt{vect} with a 3-D vector and \ttt{func} with `=' for translate to the given posiition, with `-' for to subtract the value of \ttt{vect} from the current position, and with `+' to add the value of \ttt{vect} to the current position. +Translates an entire filament in space without rescaling or rotation. Enter \ttt{vect} with a 3D vector and \ttt{func} with `=' for translate to the given posiition, with `-' for to subtract the value of \ttt{vect} from the current position, and with `+' to add the value of \ttt{vect} to the current position. For 2D simulations, \ttt{vect[2]} needs to equal 0. \item[\ttt{void}] \ttt{filRotateVertex(filamentptr fil, int seg, const double *angle, char endchar, char func);} \hfill \\ -Not written yet. This function is supposed to rotate part of the filament about one of the filament vertices by ypr angle \ttt{angle}. The character inputs give which end moves and whether it moves the angle to the given value, whether it subtracts the given value from the current value, or whether it adds the current value. +This function rotates part of the filament about the vertex that is at the front of segment \ttt{seg} by ypr angle \ttt{angle}. Enter \ttt{endchar} as \ttt{b} to have the back of the filament move and as \ttt{f} to have the front move. Enter \ttt{func} as \ttt{=} to set the current relative rotation angles (\ttt{segment->dcm}) to the input ones, as \ttt{+} for increased rotation by the given angles and as \ttt{-} for decreased rotation by these angles. + +This function makes the requested modification to the \ttt{dcm} angles and then propagates the changes either backward or forward as needed. During propagation, new segment positions are computed from neighbor positions, rather than being rotated independently, in order to avoid round-off errors. \item[\ttt{void}] \ttt{filLengthenSegment(filamentptr fil, int seg, double length, char endchar, char func);} \hfill \\ -Not written yet. This function is supposed to modify the length of a single segment. As before, the character inputs give which end moves and whether it moves the length to the given value, whether it subtracts the given value from the current value, or whether it adds the current value. +This function modifies the length of a single segment of index \ttt{seg} (not accounting for \ttt{fil->frontseg}) in filament \ttt{fil}. The \ttt{endchar} input can be \ttt{b} or \ttt{f} and gives which end of the segment moves, along with the entire rest of that end of the filament. The \ttt{func} input can be \ttt{=}, \ttt{+}, or \ttt{-} and gives whether the segment length should be made equal to \ttt{length}, increased by \ttt{length}, or decreased by \ttt{length}. This function is designed to run fairly fast, so doesn't loop over dimensions and it doesn't call functions from elsewhere. \item[\ttt{void filReverseFilament(filamentptr fil);}] \hfill \\ -Not written yet. This will reverse the sequence of segments in a filament. +NOT WRITTEN YET. This will reverse the sequence of segments in a filament. \item[\ttt{int}] \ttt{filCopyFilament(filamentptr filto, const filamentptr filfrom);} \hfill \\ This function copies all of the values from \ttt{filfrom} to \ttt{filto} except for the name, completely overwriting any prior data in \ttt{filto} in the process. Any elements in bead, segment, and monomer lists are shifted back so that they start at index number 0. Returns 0 for success, 1 for out of memory, or 2 for illegal function inputs. +\item[\ttt{filamentptr}] +\ttt{filAddFilament(filamenttypeptr filtype, const char *filname);} +\hfill \\ +This adds a filament to the filament type \ttt{filtype}, expanding the list if needed. If there is already a filament with this name, a pointer to it is returned. If \ttt{filname} is \ttt{NULL}, this generates an automatic name for the new filament; otherwise, it names the new filament to \ttt{filname}. Returns a pointer to the new filament or \ttt{NULL} for failure, arising from failure to allocate memory. + \item[\underline{Filament type}] \item[\ttt{int}] \ttt{filtypeSetParam(filamenttypeptr filtype, const char *param, int index, double value)} \hfill \\ -Sets various filament type parameters. The parameter name is \ttt{param} and it is set to the value \ttt{value}. If this parameter has multiple inidices, enter the index to be changed in \ttt{index}, or enter \ttt{index} as -1 to set all of the indices at once to the same value. Returns 0 for success or 2 for out of bounds value. The parameter options are: \ttt{stdlen} for the standard length, \ttt{stdypr} for the standard yaw-pitch-roll angles, \ttt{klen} for the stretching force constant, \ttt{kypr} for the bending force constants, \ttt{kT} for the thermodynamic energy, \ttt{treadrate} for the treadmilling rate, \ttt{viscosity} for the medium viscosity, and \ttt{beadradius} for the bead radius. +Sets various filament type parameters. The parameter name is \ttt{param} and it is set to the value \ttt{value}. If this parameter has multiple inidices, enter the index to be changed in \ttt{index}, or enter \ttt{index} as -1 to set all of the indices at once to the same value. Returns 0 for success or 2 for out of bounds value. + +\begin{longtable}[c]{llll} +\ttt{param} & Description & \ttt{index} & Range\\ +\hline +\ttt{stdlen} & standard length & n/a & $[0,\infty)$\\ +\ttt{stdypr} & standard ypr angles & 0,1,2 & $[-\pi,\pi]$\\ +\ttt{klen} & stretching force constant & n/a & $-1,[0,\infty)$\\ +\ttt{kypr} & bending force constants & 0,1,2 & $-1,[0,\infty)$\\ +\ttt{kT} & thermodynamic temperature & n/a & $[0,\infty)$\\ +\ttt{treadrate} & treadmill rate constant & n/a & $[0,\infty)$\\ +\ttt{viscosity} & medium viscosity & n/a & $(0,\infty)$\\ +\ttt{bundle} & bundle value & n/a & $(0,\infty)$\\ +\ttt{radius} & filament radius & n/a & $(0,\infty)$\\ +\ttt{facetwist} & face twist rate & n/a & $(-\infty,\infty)$ +\end{longtable} \item[\ttt{int filtypeSetColor(filamenttypeptr filtype, const double *rgba)}] \hfill \\ @@ -4328,46 +4339,62 @@ \subsection*{Function declarations.} \item[\ttt{int filtypeSetStipple(filamenttypeptr filtype, int factor, int pattern)}] \hfill \\ -Sets the stippling factor and pattern for the filament type to the entered values. Returns 0 for success or 2 if values are out of bounds. +Sets the stippling factor and pattern for the filament type to the entered values. Enter a negative number to not set the given parameter. The factor should be a positive value and the pattern should be a value between 0 and \ttt{0xFFFF}. Returns 0 for success or 2 if values are out of bounds. \item[\ttt{int filtypeSetDrawmode(filamenttypeptr filtype, enum DrawMode dm)}] \hfill \\ -Sets the filament type drawing mode to the entered value. Returns 0 for success or 2 if the input value is unrecognized. +Sets the filament type drawing mode to the entered value. Returns 0 for success or 2 if the input value is unrecognized (i.e. equal to \ttt{DMnone}). \item[\ttt{int filtypeSetShiny(filamenttypeptr filtype, double shiny)}] \hfill \\ -Sets the shininess for the filament type to the entered value. Returns 0 for success or 2 if the value is out of bounds. +Sets the shininess for the filament type to the entered value, which should be within $[0,128]$. Returns 0 for success or 2 if the value is out of bounds. \item[\ttt{int filtypeSetDynamics(filamenttypeptr filtype, enum FilamentDynamics fd)}] \hfill \\ -Sets the dynamics for the filament type to \ttt{fd}. This also sets the \ttt{isbead} element. +Sets the dynamics for the filament type to \ttt{fd}. + +\item[\ttt{int filtypeSetBiology(filamenttypeptr filtype, enum FilamentBiology fb)}] +\hfill \\ +Sets the biology for the filament type to \ttt{fb}. NOT WRITTEN YET. This function should also pre-populate as many of the other parameters as possible, which hasn't been done yet. As part of this, I should also allow the user to specify units in the \ttt{sim} structure. + +\item[\ttt{int filtypeAddFace(filamenttypeptr filtype, const char* facename)}] +\hfill \\ +Adds a face name to the filament, enlarging the list of face names if required. Returns 0 for success and -1 for failure to allocate memory. + +\item[\ttt{filamenttypeptr filAddFilamentType(simptr sim, const char *ftname)}] +\hfill \\ +Adds a new filament type to the simulation, which is named \ttt{ftname}, returning a pointer to the new filament type. Both \ttt{sim} and \ttt{ftname} need to be defined. If there is already a filament type with this name, a pointer to it is returned. Returns a pointer to the new filament type or \ttt{NULL} for failure, meaning that memory could not be allocated. \item[\underline{Filament superstructure}] \item[\ttt{void filsetcondition(filamentssptr filss, enum StructCond cond, int upgrade);}] \hfill \\ -Local. This function sets the condition of the filament superstructure. Set \ttt{upgrade} to 1 if this is an upgrade, to 0 if this is a downgrade, or to 2 to set the condition independent of its current value. +Local. This function sets the condition of the filament superstructure. Set \ttt{upgrade} to 1 if this is an upgrade, to 0 if this is a downgrade, or to 2 to set the condition independent of its current value. If the resulting condition is lower than the \ttt{sim} condition, it downgrades the \ttt{sim} condition. \item[\ttt{int filenablefilaments(simptr sim);}] \hfill \\ -Allocates a filament superstructure, adds it to the simulation structure, and sets the filament condition to \ttt{SClists}. +Allocates a filament superstructure, adds it to the simulation structure, and sets the filament condition to \ttt{SClists}. This does not create any filament types or other items that are lower in the heirarchy. -\item[\ttt{filamentptr}] -\ttt{filAddFilament(filamenttypeptr filtype, filamentptr fil, const char *filname);} +\item[\ttt{int filupdateparams(simptr sim);}] \hfill \\ -This has several models, all of which generally add a new filament to the simulation, which is named \ttt{filname}, returning a pointer to the new filament. If a filament with this name already exists, this returns a pointer to the existing filament. +Does nothing currently. This function will compute any simulation parameters necessary from user parameters. -If both \ttt{filtype} and \ttt{fil} are \ttt{NULL}, this creates a new filament structure, with its own \ttt{filname} string, and returns it; it cannot participate in simulations because it's not part of a filament type, but it can be set up. This detached filament is typically then returned to this function, now with both \ttt{filtype} and \ttt{fil} defined, and this function copies it into the filament type, freeing the extra name string in the process. Alternatively, a tidier approach is to know the filament type when creating a new filament, in which case \ttt{filtype} is defined and \ttt{fil} is \ttt{NULL}; in this case, the function adds a new filament structure to the correct place in the simulation and returns a pointer to it. Finally, there's no point in calling this function with \ttt{fil} defined and \ttt{filtype} equal to \ttt{NULL} because there's nothing sensible that it can do, so it just returns \ttt{fil}. +\item[\ttt{int filupdatelists(simptr sim);}] +\hfill \\ +Does nothing currently. This function will compute any simulation list changes necessary from user parameters. -\item[\ttt{filamenttypeptr filAddFilamentType(simptr sim, const char *ftname)}] +\item[\ttt{int filupdate(simptr sim);}] \hfill \\ -Adds a new filament type to the simulation, which is named \ttt{ftname}, returning a pointer to the new filament type. If a filament type with this name already exists, this returns a pointer to the existing filament type. +Updates filaments as needed and upgrades the condition. Returns 0 for success or other error codes passed through from \ttt{filupdatelists} and \ttt{filupdateparams}. -\item[\ttt{int}] -\ttt{filAddRandomSegments(filamentptr fil, int number, const char *xstr, const char *ystr, const char *zstr, double thickness);} + +\item[\underline{User input}] + +\item[\ttt{filamenttypeptr}] +\ttt{filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr filtype,const char *word,char *line2);} \hfill \\ -Adds \ttt{number} segments to the existing filament in \ttt{fil}. The \ttt{xstr}, \ttt{ystr}, and \ttt{zstr} entries are only considered if the filament previously had no segments; in this case, the filament starting point is at the coordinates given by them. Each entry can be a `u' to indicate uniform distributed random location in the simulation volume, or a value for the actual coordinate. Enter \ttt{thickness} with the segment thickness. +Reads one line of user input within a filament type block. \item[\ttt{filamentptr}] \ttt{filreadstring(simptr sim, ParseFilePtr pfp, filamentptr fil, const char *word, char *line2);} @@ -4378,32 +4405,29 @@ \subsection*{Function declarations.} \hfill \\ Reads multiple lines of user input from within a filament block, calling \ttt{filreadstring} with each line. -\item[\ttt{int filupdateparams(simptr sim);}] -\hfill \\ -Does nothing currently. This function will take compute any simulation parameters necessary from user parameters. -\item[\ttt{int filupdatelists(simptr sim);}] -\hfill \\ -Does nothing currently. This function will take compute any simulation list changes necessary from user parameters. +\item[\underline{Core simulation functions}] -\item[\ttt{int filsupdate(simptr sim);}] +\item[\ttt{int filSegmentXSurface(simptr sim,segmentptr segment,panelptr *pnlptr);}] \hfill \\ +Tests if segment \ttt{segment} crosses any surface in the simulation. Returns 0 if not and 1 if so. If so, then this returns a pointer to the panel that is crossed in \ttt{pnlptr}. Many more details about the crossing could be returned, such as where along the segment the crossing is, if there are multiple crossings with curved panels, etc., but these are not being returned at present. This function scans through the possible surface panels and checks each with \ttt{lineXpanel}. -\item[\underline{Core simulation functions}] - -\item[\ttt{int filMonomerXSurface(simptr sim, filamentptr fil, char endchar);}] +\item[\ttt{int filSegmentXFilament(simptr sim,segmentptr segment,filamentptr *filptr);}] \hfill \\ +Tests if segment \ttt{segment} crosses any filament in the simulation. Returns 0 if not and 1 if so. If so, then this returns a pointer to the filament that is crossed in \ttt{filptr}. More details about the crossing could be returned, such as which segment within the filament, but this is not supported at present. This function scans through all filaments in the simulation and tests for crossing for each one, meaning that the segment radius overlaps the radius of the other filament (this uses the \ttt{thk} element as the radius). This ignores possible overlaps with the self-segment, along with its nearest neighbor on either side. \item[\ttt{int}] -\ttt{filMonomerXFilament(simptr sim, filamentptr fil, char endchar, filamentptr *filptr);} +\ttt{filAddOneRandomSegment(simptr sim,filamentptr fil,const double *x,double thickness,char endchar,int constraints);} \hfill \\ +Adds one segment to the \ttt{endchar} end of \ttt{fil} using a random length and bending angle. The \ttt{x} input is ignored unless the filament has 0 current segments. The new segment has thickness \ttt{thickness}. The \ttt{constraints} input is supposed to be a collection of flags for which constraints the addition needs to account for. For now, the only \ttt{constraints} options are 0 for no constraints at all, and 1 for the new segment should not cross any surface panel. If the addition is constrained, this function tries up to \ttt{FILMAXTRIES} times to add a segment that does not disobey the constraints. Returns 0 for success, 1 for out of memory, and 2 for failure to place a segment that obeys the constraints. -\item[\ttt{void filTreadmill(simptr sim, filamentptr fil, int steps);}] +\item[\ttt{void filTreadmill(simptr sim, filamentptr fil, char endchar);}] \hfill \\ +This performs one step of treadmilling for filament \ttt{fil} at the \ttt{endchar} end. It adds one segment at this end and removes one segment from the opposite end. \item[\ttt{int filDynamics(simptr sim);}] \hfill \\ - +This will run all dynamics for filaments. It's not written yet. \end{description} diff --git a/docs/Smoldyn/SmoldynManual.pdf b/docs/Smoldyn/SmoldynManual.pdf index 495c2160..33ce885b 100644 Binary files a/docs/Smoldyn/SmoldynManual.pdf and b/docs/Smoldyn/SmoldynManual.pdf differ diff --git a/docs/Smoldyn/SmoldynManual.tex b/docs/Smoldyn/SmoldynManual.tex index d1751b1a..d1a14f3c 100644 --- a/docs/Smoldyn/SmoldynManual.tex +++ b/docs/Smoldyn/SmoldynManual.tex @@ -3530,9 +3530,13 @@ \section{Statements about filaments} Start of filament type definition block. The filament type name may be given with $name$, or it may be given afterward with the \ttt{name} statement. If the name has not been used yet for a filament type, then a new filament type is started. Between this instruction and \ttt{end\_filament\_type}, all lines need to pertain to filament types. Parameters of one filament type can be listed in multiple blocks, or parameters for many filament types can be listed in one block. +\item{\ttt{name} $name$} + +From within a filament type definition block, this switches to a new filament type and creates it if needed. + \item{\ttt{dynamics} $dynamics$} -Sets the dynamics for the filament type. Options are: ``RigidBeads'', ``RigidSegments'', ``Rouse'', ``Alberts'', ``Nedelec''. These are case-insensitive. +Sets the dynamics for the filament type. Current options are: ``none'', ``Newton''. These are case-insensitive. \item{\ttt{biology} $biology$} @@ -3570,17 +3574,13 @@ \section{Statements about filaments} Viscosity of the surrounding medium. -\item{\ttt{bead\_radius} $value$} - -Radius of the beads in these filaments. Used for hydrodynamic calculations. - \item{\ttt{standard\_length} $length$} Relaxed length of a filament segment. It can change through stretching or compression. \item{\ttt{standard\_angle} $yaw$ $pitch$ $roll$} -Relaxed angles between adjacent filament segments. When facing toward the filament's front end, yaw represents left-right bending, pitch represents up-down bending, and roll represents rotation about the filament axis. +Relaxed angles between adjacent filament segments. When facing toward the filament's front end, yaw represents left-right bending, pitch represents up-down bending, and roll represents rotation about the filament axis. For 2D simulations, only enter a single bending angle. \item{\ttt{force\_length} $value$} @@ -3600,21 +3600,21 @@ \section{Statements about filaments} \begin{description} -\item{\ttt{random\_filament} $name$ $type$ $segments$ $[x$ $y$ $z]$} +\item{\ttt{random\_filament} $name$ $type$ $segments$ $[x$ $y$ $z$ $\theta$ $\phi$ $\chi]$ $thickness$} -Create a new filament with random segments. It is named $name$, is of type $type$, has $segments$ number of segments, and, optionally, has its starting location at $(x,y,z)$. +Create a new filament with random segments. It is named $name$, is of type $type$, has $segments$ number of segments, and, optionally, has its starting location at $(x,y,z)$ and initial yaw-pitch-roll angle $(\theta, \phi, \chi)$. It also has optional thickness $thickness.$. -\item{\ttt{start\_filament} $name$} +\item{\ttt{start\_filament} $type$ $name$} -Start of filament definition block. The filament name may be given with $name$, or it may be given afterward with the \ttt{name} statement. If the name has not been used yet for a filament, then a new filament is started. Between this instruction and \ttt{end\_filament}, all lines need to pertain to filaments. Parameters of one filament can be listed in multiple blocks, or parameters for many filaments can be listed in one block. +Start of filament definition block, for a filament of type $type$ and named $name$. If the name has not been used yet for a filament of this type, then a new filament is started. Between this instruction and \ttt{end\_filament}, all lines need to pertain to filaments. Parameters of one filament can be listed in multiple blocks, or parameters for many filaments can be listed in one block. -\item{\ttt{* name} $name$} +\item{\ttt{new\_filament} $type$ $name$} -Name of the filament for editing. This statement is not required because the filament name can also be given with \ttt{start\_filament}. This statement gives the name of the current filament for editing, and creates a new filament if needed. +Defines a new filament of type $type$ that is called $name$, but does not start a filament block. This statement is largely redundant with \ttt{start\_filament}. -\item{\ttt{type} $type$} +\item{\ttt{* name} $type$ $name$} -Filament type, which should already be defined in a filament type block. +Name of the filament for editing. This statement is not required because the filament name can also be given with \ttt{start\_filament}. This statement gives the name of the current filament for editing, and creates a new filament if needed. \item{\ttt{first\_segment} $x$ $y$ $z$ $length$ $angle_0$ $angle_1$ $angle_2$ $[thickness]$} diff --git a/examples/S13_filaments/filament.txt b/examples/S13_filaments/filament.txt index 7e6df856..be7e23a2 100644 --- a/examples/S13_filaments/filament.txt +++ b/examples/S13_filaments/filament.txt @@ -1,8 +1,8 @@ # Filaments -#graphics opengl -random_seed 2 -graphics opengl_better +graphics opengl +#random_seed 2 +#graphics opengl_better dim 3 boundaries x 0 100 r @@ -13,7 +13,7 @@ species red difc red 3 -color red white +color red red time_start 0 time_stop 10000 @@ -21,7 +21,7 @@ time_step 0.01 frame_thickness 0 -mol 1 red u u u +mol 2 red u u u start_filament_type black @@ -32,18 +32,22 @@ kT 1 #treadmill_rate 1 standard_length 5 standard_angle 0 0 0 -force_length 0 +force_length 1 force_angle 100 0 0 end_filament_type -start_filament template -type black +start_filament black template random_segments 5 10 50 1 u u u +#copy fil2 +#name black fil2 +#translate + 15 15 0 end_filament -#filament template copy fil2 +filament template copy fil2 +filament fil2 translate + 15 15 0 + /* filament template copy fil3 diff --git a/examples/S13_filaments/filament1.txt b/examples/S13_filaments/filament1.txt index cdfc80f1..2b2b7da1 100644 --- a/examples/S13_filaments/filament1.txt +++ b/examples/S13_filaments/filament1.txt @@ -26,14 +26,18 @@ color black thickness 2 polygon edge kT 1 -#treadmill_rate 1 -standard_length 5 +treadmill_rate 2 +standard_length 2 standard_angle 0 0 0 -force_length 0 -force_angle 100 0 0 +force_length 10 +force_angle 5 -1 -1 end_filament_type -random_filament fil1 tread 5 10 50 1 u u u +random_filament fil1 tread 10 10 50 1 0 0 0 +random_filament fil2 tread 10 20 50 1 1 0 0 +random_filament fil3 tread 10 10 40 1 -1 0 0 +random_filament fil4 tread 10 30 50 1 2 0 0 +random_filament fil5 tread 10 30 30 1 -1 0 0 start_surface bounds diff --git a/source/Smoldyn/smoldyn.h b/source/Smoldyn/smoldyn.h index 2f2a38b5..521d5a43 100644 --- a/source/Smoldyn/smoldyn.h +++ b/source/Smoldyn/smoldyn.h @@ -756,15 +756,6 @@ typedef struct latticesuperstruct /********************************* Filaments ********************************/ -enum FilamentDynamics -{ - FDnone, - FDrigidbeads, - FDrigidsegments, - FDrouse, - FDalberts, - FDnedelec -}; enum FilamentBiology { FBactin, @@ -776,37 +767,38 @@ enum FilamentBiology FBnone }; -typedef struct beadstruct +enum FilamentDynamics { - double xyz[3]; // bead coordinates - double xyzold[3]; // bead coordinates for prior time -} * beadptr; + FDnone, + FDnewton +}; typedef struct segmentstruct { - double xyzfront[3]; // Coords. for segment front - double xyzback[3]; // Coords. for segment back - double len; // segment length - double ypr[3]; // relative ypr angles - double dcm[9]; // relative dcm - double adcm[9]; // absolute segment orientation - double thk; // thickness of segment + struct filamentstruct* fil; // owning filament + int index; // self index along filament + double xyzfront[3]; // Coords. for segment front + double xyzback[3]; // Coords. for segment back + double len; // segment length + double thk; // thickness of segment + double ypr[3]; // relative ypr angles + double dcm[9]; // relative dcm + double adcm[9]; // absolute segment orientation } * segmentptr; typedef struct filamentstruct { - struct filamenttypestruct* filtype; // filament type structure - char* filname; // filament name (reference, not owned) - int maxbs; // number of beads or segments allocated - int nbs; // number of beads or segments - int frontbs; // index of front bead or segment - beadptr* beads; // array of beads if any + struct filamenttypestruct* filtype; // owning filament type + char* filname; // filament name (ref, not owned) + int maxseg; // number of segments allocated + int nseg; // number of segments + int frontseg; // index of front segment segmentptr* segments; // array of segments if any - struct filamentstruct* frontend; // what front attaches to if anything - struct filamentstruct* backend; // what back attaches to if anything - int maxbranch; // max number of branches off this filament - int nbranch; // number of branches off this filament - int* branchspots; // list of bead or segments where branches are + struct filamentstruct* frontend; // what front attaches to + struct filamentstruct* backend; // what back attaches to + int maxbranch; // max branches off this filament + int nbranch; // num branches off this filament + int* branchspots; // segments where branches are struct filamentstruct** branches; // list of branching filaments int maxmonomer; // number of monomers allocated int nmonomer; // number of monomers @@ -817,11 +809,10 @@ typedef struct filamentstruct typedef struct filamenttypestruct { struct filamentsuperstruct* filss; // filament superstructure - char* ftname; // filament type name (reference, not owned) + char* ftname; // fil type name (ref, not owned) + enum FilamentBiology biology; // biological name for filaments enum FilamentDynamics dynamics; // dynamics for the filament - int isbead; // 1 for bead model, 0 for segments - enum FilamentBiology biology; // Biological name for filament type - double bundlevalue; // number of microfilaments in bundle + double bundlevalue; // microfilaments in bundle double color[4]; // filament color double edgepts; // thickness of edge for drawing unsigned int edgestipple[2]; // edge stippling [factor, pattern] @@ -831,16 +822,17 @@ typedef struct filamenttypestruct double stdypr[3]; // minimum energy bend angle double klen; // force constant for length double kypr[3]; // force constant for angle - double kT; // thermodynamic temperature, [0,inf) + double kT; // thermodynamic temp, [0,inf) double treadrate; // treadmilling rate constant double viscosity; // viscosity - double filradius; // bead or segment radius - int maxface; // number of filament faces allocated + double filradius; // segment radius + int maxface; // filament faces allocated int nface; // number of filament faces char** facename; // list of face names - double facetwist; // twisting rate of faces along filament + double facetwist; // twisting rate of faces int maxfil; // maximum number of filaments int nfil; // actual number of filaments + int autonamenum; // next number for autoname filamentptr* fillist; // list of filaments char** filnames; // names of filaments } * filamenttypeptr; @@ -855,44 +847,6 @@ typedef struct filamentsuperstruct filamenttypeptr* filtypes; // list of filament types } * filamentssptr; -/* OLD FILAMENTS -typedef struct filamentstruct -{ - struct filamentsuperstruct* filss; // filament superstructure - char* fname; // filament name - double color[4]; // filament color - double edgepts; // thickness of edge for drawing - unsigned int edgestipple[2]; // edge stippling [factor, pattern] - enum DrawMode drawmode; // polygon drawing mode - double shiny; // shininess - int maxseg; // number of segments allocated - int nseg; // number of segments - int front; // front index - int back; // back index - double** sxyz; // Coords. for segment ends [seg][3] - double* slen; // segment length [seg] - double** sypr; // relative ypr angles [seg][3] - double** sdcm; // relative dcm [seg][9] - double** sadcm; // absolute segment orientation [seg][9] - double* sthk; // thickness of segment [nmax], [0,inf) - double stdlen; // minimum energy segment length - double stdypr[3]; // minimum energy bend angle - double klen; // force constant for length - double kypr[3]; // force constant for angle - double kT; // thermodynamic temperature, [0,inf) - double treadrate; // treadmilling rate constant -} * filamentptr; - -typedef struct filamentsuperstruct -{ - enum StructCond condition; // structure condition - struct simstruct* sim; // simulation structure - int maxfil; // maximum number of filaments - int nfil; // actual number of filaments - char** fnames; // filament names - filamentptr* fillist; // list of filaments -} * filamentssptr; -*/ /******************************** BioNetGen *********************************/ diff --git a/source/Smoldyn/smoldynfuncs.h b/source/Smoldyn/smoldynfuncs.h index ddad13ab..79c7a2e9 100644 --- a/source/Smoldyn/smoldynfuncs.h +++ b/source/Smoldyn/smoldynfuncs.h @@ -410,17 +410,17 @@ void filssoutput(simptr sim); int filcheckparams(simptr sim,int *warnptr); // filament manipulation -int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr,const char *thtstr,const char *phistr,const char *chistr,double thickness); -int filAddRandomBeads(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr); +int filGetFilIndex(simptr sim,const char *name,int *ftptr); +filamentptr filAddFilament(filamenttypeptr filtype,const char *filname); +int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr,const char *thtstr, const char* phistr,const char* chistr,double thickness); // structure set up int filenablefilaments(simptr sim); -filamentptr filAddFilament(filamenttypeptr filtype,filamentptr fil,const char *filname); filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr filtype,const char *word,char *line2); filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamenttypeptr filtype,const char *word,char *line2); int filloadtype(simptr sim,ParseFilePtr *pfpptr,char *line2); -int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2,filamenttypeptr filtype); -int filsupdate(simptr sim); +int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2); +int filupdate(simptr sim); // core simulation functions int filDynamics(simptr sim); diff --git a/source/Smoldyn/smolfilament.c b/source/Smoldyn/smolfilament.c index 4221aec2..2ff29afd 100644 --- a/source/Smoldyn/smolfilament.c +++ b/source/Smoldyn/smolfilament.c @@ -19,6 +19,9 @@ #include "smoldynfuncs.h" +#define FILMAXTRIES 100 + + /******************************************************************************/ /********************************** Filaments *********************************/ /******************************************************************************/ @@ -29,60 +32,76 @@ /******************************************************************************/ // enumerated types - +char *filFB2string(enum FilamentBiology fb,char *string); +enum FilamentBiology filstring2FB(const char *string); +char *filFD2string(enum FilamentDynamics fd,char *string); +enum FilamentDynamics filstring2FD(const char *string); // low level utilities -double filStretchEnergy(filamentptr fil,int seg1,int seg2); -double filBendEnergy(filamentptr fil,int seg1,int seg2); +double filRandomLength(const filamenttypeptr filtype,double thickness,double sigmamult); +double *filRandomAngle(const filamenttypeptr filtype,int dim,int n,double thickness,double sigmamult,double *angle); + +// computations on filaments +double filStretchEnergy(const filamentptr fil,int seg1,int seg2); +double filBendEnergy(const filamentptr fil,int seg1,int seg2); // memory management -void beadfree(beadptr bead); +segmentptr segmentalloc(); void segmentfree(segmentptr segment); -void filamenttypefree(filamenttypeptr filtype); -filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxmonomer); +filamentptr filalloc(filamentptr fil,int maxseg,int maxbranch,int maxmonomer); void filfree(filamentptr fil); - +filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface); +void filamenttypefree(filamenttypeptr filtype); +filamentssptr filssalloc(filamentssptr filss,int maxtype); +void filssfree(filamentssptr filss); // data structure output - -// structure set up +void filoutput(const filamentptr fil); +void filtypeoutput(const filamenttypeptr filtype); +void filssoutput(const simptr sim); +void filwrite(const simptr sim,FILE *fptr); +int filcheckparams(const simptr sim,int *warnptr); + +// filaments +void filUpdateSegmentIndices(filamentptr fil); +void filArrayShift(filamentptr fil,int shift); +int filAddSegment(filamentptr fil,const double *x,double length,const double *angle,double thickness,char endchar); +int filRemoveSegment(filamentptr fil,char endchar); +void filTranslate(filamentptr fil,const double *vect,char func); +void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func); +void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func); +void filReverseFilament(filamentptr fil); +int filCopyFilament(filamentptr filto,const filamentptr filfrom); + +// filament type +int filtypeSetParam(filamenttypeptr filtype,const char *param,int index,double value); int filtypeSetColor(filamenttypeptr filtype,const double *rgba); int filtypeSetEdgePts(filamenttypeptr filtype,double value); int filtypeSetStipple(filamenttypeptr filtype,int factor,int pattern); int filtypeSetDrawmode(filamenttypeptr filtype,enum DrawMode dm); int filtypeSetShiny(filamenttypeptr filtype,double shiny); +int filtypeSetDynamics(filamenttypeptr filtype,enum FilamentDynamics fd); int filtypeSetBiology(filamenttypeptr filtype,enum FilamentBiology fb); +int filtypeAddFace(filamenttypeptr filtype,const char* facename); +filamenttypeptr filAddFilamentType(simptr sim,const char *ftname); +// filament superstructure +void filsetcondition(filamentssptr filss,enum StructCond cond,int upgrade); +int filenablefilaments(simptr sim); int filupdateparams(simptr sim); int filupdatelists(simptr sim); +int filupdate(simptr sim); - -/******************************************************************************/ -/********************************* enumerated types ***************************/ -/******************************************************************************/ +// user input -/* filFD2string */ -char *filFD2string(enum FilamentDynamics fd,char *string) { - if(fd==FDnone) strcpy(string,"none"); - else if(fd==FDrouse) strcpy(string,"rouse"); - else if(fd==FDalberts) strcpy(string,"alberts"); - else if(fd==FDnedelec) strcpy(string,"nedelec"); - else strcpy(string,"none"); - return string; } -/* filstring2FD */ -enum FilamentDynamics filstring2FD(const char *string) { - enum FilamentDynamics ans; - if(strbegin(string,"none",0)) ans=FDnone; - else if(strbegin(string,"rouse",0)) ans=FDrouse; - else if(strbegin(string,"alberts",0)) ans=FDalberts; - else if(strbegin(string,"nedelec",0)) ans=FDnedelec; - else ans=FDnone; - return ans; } +/******************************************************************************/ +/********************************* enumerated types ***************************/ +/******************************************************************************/ /* filFB2string */ @@ -92,7 +111,8 @@ char *filFB2string(enum FilamentBiology fb,char *string) { else if(fb==FBintermediate) strcpy(string,"intermediate"); else if(fb==FBdsDNA) strcpy(string,"dsDNA"); else if(fb==FBssDNA) strcpy(string,"ssDNA"); - else strcpy(string,"other"); + else if(fb==FBother) strcpy(string,"other"); + else strcpy(string,"none"); return string; } @@ -110,6 +130,22 @@ enum FilamentBiology filstring2FB(const char *string) { return ans; } +/* filFD2string */ +char *filFD2string(enum FilamentDynamics fd,char *string) { + if(fd==FDnewton) strcpy(string,"newton"); + else strcpy(string,"none"); + return string; } + + +/* filstring2FD */ +enum FilamentDynamics filstring2FD(const char *string) { + enum FilamentDynamics ans; + + if(strbegin(string,"none",0)) ans=FDnone; + else if(strbegin(string,"newton",0)) ans=FDnewton; + else ans=FDnone; + return ans; } + /******************************************************************************/ /****************************** low level utilities ***************************/ @@ -121,7 +157,7 @@ double filRandomLength(const filamenttypeptr filtype,double thickness,double sig double lsigma; double length; - if(filtype->klen<=0) return filtype->stdlen; + if(filtype->klen<0) return filtype->stdlen; lsigma=sigmamult*sqrt(filtype->kT/(thickness*filtype->klen)); length=0; while(length<=0) @@ -130,18 +166,32 @@ double filRandomLength(const filamenttypeptr filtype,double thickness,double sig /* filRandomAngle */ -double *filRandomAngle(const filamenttypeptr filtype,double thickness,double *angle,double sigmamult) { - double sigma; +double *filRandomAngle(const filamenttypeptr filtype,int dim,int n,double thickness,double sigmamult,double *angle) { int d; - for(d=0;d<3;d++){ - if(filtype->kypr[d]>0) - sigma=sigmamult*sqrt(filtype->kT/(thickness*filtype->kypr[d])); - else if(filtype->kypr[d]==0) - sigma=unirandOCD(-PI,PI); + if(n>0 && dim==3) { // 3D not first segment + for(d=0;d<3;d++){ + if(filtype->kypr[d]<0) + angle[d]=filtype->stdypr[d]; + else if(filtype->kypr[d]==0) + angle[d]=unirandOCD(-PI,PI); + else + angle[d]=filtype->stdypr[d]+sigmamult*sqrt(filtype->kT/(thickness*filtype->kypr[d]))*gaussrandD(); }} + else if(n>0) { // 2D not first segment + if(filtype->kypr[0]<0) + angle[0]=filtype->stdypr[0]; + else if(filtype->kypr[0]==0) + angle[0]=unirandOCD(-PI,PI); else - sigma=0; - angle[d]=filtype->stdypr[d]+(sigma>0?sigma*gaussrandD():0); } + angle[0]=filtype->stdypr[0]+sigmamult*sqrt(filtype->kT/(thickness*filtype->kypr[0]))*gaussrandD(); + angle[1]=angle[2]=0; } + else if(dim==3) { // 3D first segment + angle[0]=thetarandCCD(); + angle[1]=unirandCOD(0,2*PI); + angle[2]=unirandCOD(0,2*PI); } + else { // 2D first segment + angle[0]=unirandOCD(0,2*PI); + angle[1]=angle[2]=0; } return angle; } @@ -152,9 +202,9 @@ double *filRandomAngle(const filamenttypeptr filtype,double thickness,double *an // filStretchEnergy double filStretchEnergy(const filamentptr fil,int seg1,int seg2) { - double thk,klen,len,stdlen,energy,*xyz,*xyzplus; + double thk,klen,len,stdlen,energy; filamenttypeptr filtype; - int bs; + int seg; filtype=fil->filtype; stdlen=filtype->stdlen; @@ -162,49 +212,37 @@ double filStretchEnergy(const filamentptr fil,int seg1,int seg2) { if(klen<=0) return 0; energy=0; - if(seg1==-1) seg1=fil->frontbs; - if(seg2==-1) seg2=fil->frontbs+fil->nbs; - if(filtype->isbead) { - for(bs=seg1;bsbeads[bs]->xyz; - xyzplus=fil->beads[bs+1]->xyz; - len=sqrt((xyzplus[0]-xyz[0])*(xyzplus[0]-xyz[0])+(xyzplus[1]-xyz[1])*(xyzplus[1]-xyz[1])+(xyzplus[2]-xyz[2])*(xyzplus[2]-xyz[2])); - energy+=0.5*klen*(len-stdlen)*(len-stdlen); }} - else { - for(bs=seg1;bssegments[bs]->thk; - len=fil->segments[bs]->len; - energy+=0.5*thk*klen*(len-stdlen)*(len-stdlen); }} + if(seg1==-1) seg1=fil->frontseg; + if(seg2==-1) seg2=fil->frontseg+fil->nseg; + for(seg=seg1;segsegments[seg]->thk; + len=fil->segments[seg]->len; + energy+=0.5*thk*klen*(len-stdlen)*(len-stdlen); } return energy; } // filBendEnergy double filBendEnergy(const filamentptr fil,int seg1,int seg2) { - double *kypr,*stdypr,*ypr,thk,thkminus,energy; + double *kypr,*stdypr,*ypr,thk,energy; filamenttypeptr filtype; segmentptr segptr,segptrminus; - double tk; int seg; filtype=fil->filtype; - if(filtype->isbead) return 0; - kypr=filtype->kypr; stdypr=filtype->stdypr; energy=0; - if(seg1==-1) seg1=fil->frontbs; - if(seg2==-1) seg2=fil->frontbs+fil->nbs; - for(seg=seg1+1;segfrontseg; + if(seg2==-1) seg2=fil->frontseg+fil->nseg-1; + for(seg=seg1+1;seg<=seg2;seg++) { segptr = fil->segments[seg]; segptrminus = fil->segments[seg-1]; ypr = segptr->ypr; - thk = segptr->thk; - thkminus = segptrminus->thk; - tk=0.5*(thkminus+thk); - if(kypr[0]>0) energy+=0.5*tk*kypr[0]*(ypr[0]-stdypr[0])*(ypr[0]-stdypr[0]); - if(kypr[1]>0) energy+=0.5*tk*kypr[1]*(ypr[1]-stdypr[1])*(ypr[1]-stdypr[1]); - if(kypr[2]>0) energy+=0.5*tk*kypr[2]*(ypr[2]-stdypr[2])*(ypr[2]-stdypr[2]); } + thk=0.5*(segptrminus->thk+segptr->thk); + if(kypr[0]>0) energy+=0.5*thk*kypr[0]*(ypr[0]-stdypr[0])*(ypr[0]-stdypr[0]); + if(kypr[1]>0) energy+=0.5*thk*kypr[1]*(ypr[1]-stdypr[1])*(ypr[1]-stdypr[1]); + if(kypr[2]>0) energy+=0.5*thk*kypr[2]*(ypr[2]-stdypr[2])*(ypr[2]-stdypr[2]); } return energy; } @@ -212,38 +250,20 @@ double filBendEnergy(const filamentptr fil,int seg1,int seg2) { /******************************* memory management ****************************/ /******************************************************************************/ -/* beadalloc */ -beadptr beadalloc() { - beadptr bead; - - CHECKMEM(bead=(beadptr) malloc(sizeof(struct beadstruct))); - bead->xyz[0]=bead->xyz[1]=bead->xyz[2]=0; - bead->xyzold[0]=bead->xyzold[1]=bead->xyzold[2]=0; - return bead; - - failure: - return NULL; } - - -/* beadfree */ -void beadfree(beadptr bead) { - if(!bead) return; - free(bead); - return; } - /* segmentalloc */ segmentptr segmentalloc() { segmentptr segment; CHECKMEM(segment=(segmentptr) malloc(sizeof(struct segmentstruct))); + segment->index=0; segment->xyzfront[0]=segment->xyzfront[1]=segment->xyzfront[2]=0; segment->xyzback[0]=segment->xyzback[1]=segment->xyzback[2]=0; segment->len=0; + segment->thk=1; segment->ypr[0]=segment->ypr[1]=segment->ypr[2]=0; Sph_One2Dcm(segment->dcm); Sph_One2Dcm(segment->adcm); - segment->thk=1; return segment; failure: @@ -258,9 +278,8 @@ void segmentfree(segmentptr segment) { /* filalloc */ -filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxmonomer) { - int bs,br,mn; - beadptr *newbeads; +filamentptr filalloc(filamentptr fil,int maxseg,int maxbranch,int maxmonomer) { + int seg,br,mn; segmentptr *newsegments; int *newbranchspots; filamentptr *newbranches; @@ -270,10 +289,9 @@ filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxm CHECKMEM(fil=(filamentptr) malloc(sizeof(struct filamentstruct))); fil->filtype=NULL; fil->filname=NULL; - fil->maxbs=0; - fil->nbs=0; - fil->frontbs=0; - fil->beads=NULL; + fil->maxseg=0; + fil->nseg=0; + fil->frontseg=0; fil->segments=NULL; fil->frontend=NULL; fil->backend=NULL; @@ -286,28 +304,18 @@ filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxm fil->frontmonomer=0; fil->monomers=NULL; } - if(maxbs>fil->maxbs) { - if(isbead) { - CHECKMEM(newbeads=(beadptr*) calloc(maxbs,sizeof(struct beadstruct))); - for(bs=0;bsmaxbs;bs++) - newbeads[bs]=fil->beads[bs]; - for(;bsbeads); - fil->beads=newbeads; } - else { - CHECKMEM(newsegments=(segmentptr*) calloc(maxbs,sizeof(struct segmentstruct))); - for(bs=0;bsmaxbs;bs++) - newsegments[bs]=fil->segments[bs]; - for(;bssegments); - fil->segments=newsegments; } - fil->maxbs=maxbs; } + if(maxseg>fil->maxseg) { + CHECKMEM(newsegments=(segmentptr*) calloc(maxseg,sizeof(struct segmentstruct))); + for(seg=0;segmaxseg;seg++) + newsegments[seg]=fil->segments[seg]; + for(;segsegments); + fil->segments=newsegments; + fil->maxseg=maxseg; + filUpdateSegmentIndices(fil); } if(maxbranch>fil->maxbranch) { CHECKMEM(newbranchspots=(int*) calloc(maxbranch,sizeof(int))); @@ -340,16 +348,12 @@ filamentptr filalloc(filamentptr fil,int isbead,int maxbs,int maxbranch,int maxm /* filfree */ void filfree(filamentptr fil) { - int bs; + int seg; if(!fil) return; - if(fil->beads) { - for(bs=0;bsmaxbs;bs++) - beadfree(fil->beads[bs]); - free(fil->beads); } if(fil->segments) { - for(bs=0;bsmaxbs;bs++) - segmentfree(fil->segments[bs]); + for(seg=0;segmaxseg;seg++) + segmentfree(fil->segments[seg]); free(fil->segments); } free(fil->branchspots); @@ -369,9 +373,8 @@ filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface CHECKMEM(filtype=(filamenttypeptr) malloc(sizeof(struct filamenttypestruct))); filtype->filss=NULL; filtype->ftname=NULL; + filtype->biology=FBnone; filtype->dynamics=FDnone; - filtype->isbead=0; - filtype->biology=FBother; filtype->bundlevalue=1; filtype->color[0]=filtype->color[1]=filtype->color[2]=0; @@ -386,7 +389,7 @@ filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface filtype->stdypr[0]=filtype->stdypr[1]=filtype->stdypr[2]=0; filtype->klen=1; filtype->kypr[0]=filtype->kypr[1]=filtype->kypr[2]=1; - filtype->kT=0; + filtype->kT=1; filtype->treadrate=0; filtype->viscosity=1; filtype->filradius=1; @@ -398,6 +401,7 @@ filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface filtype->maxfil=0; filtype->nfil=0; + filtype->autonamenum=0; filtype->fillist=NULL; filtype->filnames=NULL; } @@ -411,7 +415,7 @@ filamenttypeptr filamenttypealloc(filamenttypeptr filtype,int maxfil,int maxface newfillist[f]=filtype->fillist[f]; newfilnames[f]=filtype->filnames[f]; } for(;ffiltype=filtype; newfillist[f]->filname=newfilnames[f]; } @@ -528,16 +532,14 @@ void filssfree(filamentssptr filss) { /* filoutput */ void filoutput(const filamentptr fil) { - int bs,br,mn,isbead,dim; + int seg,br,mn,dim; simptr sim; - beadptr bead; segmentptr segment; if(!fil) { simLog(NULL,2," NULL filament\n"); return; } - isbead=fil->filtype?fil->filtype->isbead:0; if(fil->filtype && fil->filtype->filss) { sim=fil->filtype->filss->sim; dim=sim->dim; } @@ -546,40 +548,32 @@ void filoutput(const filamentptr fil) { dim=3; } simLog(sim,2," Filament: %s\n",fil->filname); - simLog(sim,1," allocated beads or segments: %i\n",fil->maxbs); - simLog(sim,2," number of %s: %i\n",isbead?"beads":"segments",fil->nbs); - simLog(sim,1," front index: %i\n",fil->frontbs); - - if(isbead) - simLog(sim,2," bead, position\n"); - else - simLog(sim,2," segment, thickness, length, front position, relative angle\n"); - for(bs=0;bsnbs;bs++) { - if(isbead) { - bead=fil->beads[bs+fil->frontbs]; - if(dim==2) - simLog(sim,bs>5?1:2," %i pos.=(%1.3g %1.3g)\n",bs,bead->xyz[0],bead->xyz[1]); - else - simLog(sim,bs>5?1:2," %i x=(%1.3g %1.3g %1.3g)\n",bs,bead->xyz[0],bead->xyz[1],bead->xyz[2]); } - else { - segment=fil->segments[bs+fil->frontbs]; - if(dim==2) - simLog(sim,bs>5?1:2," %i thick=%1.3g, length=%1.3g, front pos.=(%1.3g %1.3g), rel. angle=%1.3g\n",bs,segment->thk,segment->len,segment->xyzfront[0],segment->xyzfront[1],segment->ypr[0]); - else - simLog(sim,bs>5?1:2," %i thick=%1.3g, length=%1.3g, front pos.=(%1.3g %1.3g %1.3g), rel. angle=(%1.3g %1.3g %1.3g)\n",bs,segment->thk,segment->len,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); }} + simLog(sim,1," type: %s\n",fil->filtype?fil->filtype->ftname:"None (assuming dim=3)"); + simLog(sim,1," allocated segments: %i\n",fil->maxseg); + simLog(sim,2," number of segments: %i\n",fil->nseg); + simLog(sim,1," front index: %i\n",fil->frontseg); + + simLog(sim,2," segment, length, thickness, front position, relative angle:\n"); + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg+fil->frontseg]; + if(dim==2) + simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g), rel. angle=%1.3g\n",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->ypr[0]); + else + simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g %1.3g), rel. angle=(%1.3g %1.3g %1.3g)\n",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); } if(fil->frontend) simLog(sim,2," front branched from: %s\n",fil->frontend->filname); if(fil->backend) simLog(sim,2," back branched from: %s\n",fil->backend->filname); - if(fil->nbranch>0) { - simLog(sim,2," number of branches: %i\n",fil->nbranch); - for(br=0;brnbranch;br++) - simLog(sim,2," %s at %i\n",fil->branches[br]->filname,fil->branchspots[br]); } + simLog(sim,1," allocated branches: %i\n",fil->maxbranch); + simLog(sim,fil->nbranch>0?2:1," number of branches: %i\n",fil->nbranch); + for(br=0;brnbranch;br++) + simLog(sim,2," %s at %i\n",fil->branches[br]->filname,fil->branchspots[br]); + + simLog(sim,1," monomer codes: %i of %i allocated,",fil->nmonomer,fil->maxmonomer); + simLog(sim,1," starting index: %i\n",fil->frontmonomer); if(fil->nmonomer) { - simLog(sim,1," monomer codes: %i of %i allocated,",fil->nmonomer,fil->maxmonomer); - simLog(sim,1," starting index: %i\n",fil->frontmonomer); simLog(sim,2," monomer code: "); for(mn=0;mnnmonomer;mn++) simLog(sim,2,"%c",fil->monomers[mn]); @@ -596,7 +590,7 @@ void filoutput(const filamentptr fil) { /* filtypeoutput */ void filtypeoutput(const filamenttypeptr filtype) { char string[STRCHAR]; - int isbead,fc,f,dim; + int fc,f,dim; simptr sim; if(!filtype) { @@ -609,11 +603,11 @@ void filtypeoutput(const filamenttypeptr filtype) { else { sim=NULL; dim=3; } - isbead=filtype->isbead; simLog(sim,2," Filament type: %s\n",filtype->ftname); - simLog(sim,2," dynamics: %s using %s\n",filFD2string(filtype->dynamics,string),isbead?"beads":"segments"); + simLog(sim,1," superstructure: %s\n",filtype->filss?"assigned":"missing (assuming dim=3)"); simLog(sim,2," biology: %s\n",filFB2string(filtype->biology,string)); + simLog(sim,2," dynamics: %s\n",filFD2string(filtype->dynamics,string)); simLog(sim,filtype->bundlevalue!=1?2:1," bundle value: %g\n",filtype->bundlevalue); simLog(sim,2," color: %g %g %g %g\n",filtype->color[0],filtype->color[1],filtype->color[2],filtype->color[3]); @@ -621,26 +615,25 @@ void filtypeoutput(const filamenttypeptr filtype) { if(filtype->edgestipple[1]!=0xFFFF) simLog(sim,2," edge stippling: %ui %X\n",filtype->edgestipple[0],filtype->edgestipple[1]); if(filtype->shiny!=0) simLog(sim,2," shininess: %g\n",filtype->shiny); - simLog(sim,2," %s length: %g\n",filtype->klen>0?"standard":"fixed",filtype->stdlen); + simLog(sim,2," %s length: %g\n",filtype->klen>=0?"standard":"fixed",filtype->stdlen); if(filtype->klen>0) simLog(sim,2," length force constant: %g\n",filtype->klen); - if(!isbead) { - if(dim==2) { - simLog(sim,2," standard angle: %g\n",filtype->stdypr[0]); - simLog(sim,2," bending force constant: %g\n",filtype->kypr[0]); } - else { - simLog(sim,2," standard angles: %g, %g, %g\n",filtype->stdypr[0],filtype->stdypr[1],filtype->stdypr[2]); - simLog(sim,2," bending force constants: %g, %g, %g\n",filtype->kypr[0],filtype->kypr[1],filtype->kypr[2]); }} + if(dim==2) { + simLog(sim,2," %s angle: %g\n",filtype->kypr[0]>=0?"standard":"fixed",filtype->stdypr[0]); + simLog(sim,2," bending force constant: %g\n",filtype->kypr[0]); } + else { + simLog(sim,2," standard angles: %g, %g, %g\n",filtype->stdypr[0],filtype->stdypr[1],filtype->stdypr[2]); + simLog(sim,2," bending force constants: %g, %g, %g\n",filtype->kypr[0],filtype->kypr[1],filtype->kypr[2]); } simLog(sim,2," kT: %g\n",filtype->kT); simLog(sim,filtype->treadrate>0?2:1," treadmilling rate: %g\n",filtype->treadrate); simLog(sim,2," viscosity: %g\n",filtype->viscosity); - simLog(sim,2," %s radius: %g\n",isbead?"bead":"segment",filtype->filradius); + simLog(sim,2," filament radius: %g\n",filtype->filradius); if(filtype->nface>0) { simLog(sim,2," %i faces with twist of %g:",filtype->nface,filtype->facetwist); for(fc=0;fcnface;fc++) - simLog(sim,2," %s",filtype->facename[fc]); + simLog(sim,2," %s,",filtype->facename[fc]); simLog(sim,2,"\n"); } simLog(sim,2," %i filaments:\n",filtype->nfil); @@ -705,168 +698,115 @@ int filcheckparams(const simptr sim,int *warnptr) { /******************************************************************************/ -/****************************** filament manipulation *************************/ +/********************************** Filaments *****************************/ /******************************************************************************/ +/* filUpdateSegmentIndices */ +void filUpdateSegmentIndices(filamentptr fil) { + int i; + + for(i=0;inseg;i++) + fil->segments[i+fil->frontseg]->index=i; + return; } + + +/* filGetFilIndex */ +int filGetFilIndex(simptr sim,const char *name,int *ftptr) { + int f,f1,ft,ft1; + filamenttypeptr filtype; + + f1=-1; + ft1=-1; + for(ft=0;ftfilss->ntype;ft++) { + filtype=sim->filss->filtypes[ft]; + f=stringfind(filtype->filnames,filtype->nfil,name); + if(f>=0 && f1>=0) return -2; + else if(f>=0) { + f1=f; + ft1=ft; }} + if(f1>=0) *ftptr=ft1; + return f1; } + + /* filArrayShift */ void filArrayShift(filamentptr fil,int shift) { - int i,frontbs,backbs; - int isbead; - beadptr nullbead; + int i,frontseg,backseg; segmentptr nullsegment; - isbead=fil->filtype->isbead; if(!shift) - shift=(fil->maxbs-fil->nbs)/2-fil->frontbs; + shift=(fil->maxseg-fil->nseg)/2-fil->frontseg; - frontbs=fil->frontbs; - backbs=fil->frontbs+fil->nbs; + frontseg=fil->frontseg; + backseg=fil->frontseg+fil->nseg; if(shift>0) { - if(backbs+shift>fil->maxbs) shift=fil->maxbs-backbs; - if(isbead) { - for(i=backbs+shift-1;i>=frontbs+shift;i--) { // i is destination - nullbead = fil->beads[i]; - fil->beads[i] = fil->beads[i-shift]; - fil->beads[i-shift] = nullbead; }} - else { - for(i=backbs+shift-1;i>=frontbs+shift;i--) { - nullsegment = fil->segments[i]; - fil->segments[i]=fil->segments[i-shift]; - fil->segments[i-shift] = nullsegment; }} - fil->frontbs+=shift; } + if(backseg+shift>fil->maxseg) shift=fil->maxseg-backseg; + for(i=backseg+shift-1;i>=frontseg+shift;i--) { + nullsegment = fil->segments[i]; + fil->segments[i]=fil->segments[i-shift]; + fil->segments[i-shift] = nullsegment; } + fil->frontseg+=shift; } else if(shift<0) { shift=-shift; // now shift is positive - if(frontbs-shift<0) shift=frontbs; - if(isbead) { - for(i=frontbs-shift;ibeads[i]; - fil->beads[i]=fil->beads[i+shift]; - fil->beads[i+shift] = nullbead;}} - else { - for(i=frontbs-shift;isegments[i]; - fil->segments[i]=fil->segments[i+shift]; - fil->segments[i+shift] = nullsegment; }} - fil->frontbs-=shift; } + if(frontseg-shift<0) shift=frontseg; + for(i=frontseg-shift;isegments[i]; + fil->segments[i]=fil->segments[i+shift]; + fil->segments[i+shift] = nullsegment; } + fil->frontseg-=shift; } return; } -/* filAddBead */ -int filAddBead(filamentptr fil,const double *x,double length,char endchar) { - int seg,dim; - double xyz[3]; - filamenttypeptr filtype; - beadptr bead,beadminus,beadplus; - - filtype=fil->filtype; - dim=filtype->filss->sim->dim; - if(!filtype->isbead) return -2; // can't add a bead to a non-bead filament type - - if(fil->nbs==fil->maxbs) { - fil=filalloc(fil,1,fil->maxbs*2+1,0,0); - if(!fil) return -1; } // out of memory - - if(endchar=='b') { // add to back end - if(fil->frontbs+fil->nbs==fil->maxbs) - filArrayShift(fil,0); - seg=fil->frontbs+fil->nbs; - bead=fil->beads[seg]; - if(fil->nbs==0) { - bead->xyz[0]=bead->xyzold[0]=x[0]; - bead->xyz[1]=bead->xyzold[1]=x[1]; - bead->xyz[2]=bead->xyzold[2]=x[2]; } // xyzold = xyz = x - else { - beadminus=fil->beads[seg-1]; - if(dim==2) { - circlerandD(xyz,length); - xyz[2]=0; } - else - sphererandCCD(xyz,length,length); - bead->xyz[0]=bead->xyzold[0]=beadminus->xyz[0]+xyz[0]; - bead->xyz[1]=bead->xyzold[1]=beadminus->xyz[1]+xyz[1]; - bead->xyz[2]=bead->xyzold[2]=beadminus->xyz[2]+xyz[2]; } - fil->nbs++; } - - else { // add to front end - if(fil->frontbs==0) filArrayShift(fil,0); - if(fil->frontbs==0) filArrayShift(fil,1); // used if nbs=maxbs-1 - seg=fil->frontbs-1; - bead=fil->beads[seg]; - if(fil->nbs==0) { - bead->xyz[0]=bead->xyzold[0]=x[0]; - bead->xyz[1]=bead->xyzold[1]=x[1]; - bead->xyz[2]=bead->xyzold[2]=x[2]; } - else { - beadplus=fil->beads[seg+1]; - if(dim==2) { - circlerandD(xyz,length); - xyz[2]=0; } - else - sphererandCCD(xyz,length,length); - bead->xyz[0]=bead->xyzold[0]=beadplus->xyz[0]+xyz[0]; - bead->xyz[1]=bead->xyzold[1]=beadplus->xyz[1]+xyz[1]; - bead->xyz[2]=bead->xyzold[2]=beadplus->xyz[2]+xyz[2]; } - fil->frontbs--; - fil->nbs++; } - - return 0; } - - /* filAddSegment */ int filAddSegment(filamentptr fil,const double *x,double length,const double *angle,double thickness,char endchar) { - int seg,dim; - UNUSED(dim); - filamenttypeptr filtype; + int seg; segmentptr segment,segmentminus,segmentplus; - filtype=fil->filtype; - dim=filtype->filss->sim->dim; - if(filtype->isbead) return -2; // can't add a segment to a bead filament type - - if(fil->nbs==fil->maxbs) { - fil=filalloc(fil,0,fil->maxbs*2+1,0,0); - if(!fil) return -1; } // out of memory + if(fil->nseg==fil->maxseg) { + fil=filalloc(fil,fil->maxseg*2+1,0,0); + if(!fil) return 1; } // out of memory if(endchar=='b') { - if(fil->frontbs+fil->nbs==fil->maxbs) filArrayShift(fil,0); - seg=fil->frontbs+fil->nbs; - segment=fil->segments[seg]; + if(fil->frontseg+fil->nseg==fil->maxseg) filArrayShift(fil,0); + seg=fil->nseg; + segment=fil->segments[seg+fil->frontseg]; segment->len=length; segment->thk=thickness; + segment->index=seg; Sph_Xyz2Xyz(angle,segment->ypr); // ypr = angle Sph_Xyz2Dcm(angle,segment->dcm); // A = Dcm(angle) - if(fil->nbs==0) { + if(fil->nseg==0) { segment->xyzfront[0]=x[0]; // x_i = input value segment->xyzfront[1]=x[1]; segment->xyzfront[2]=x[2]; Sph_Dcm2Dcm(segment->dcm,segment->adcm); } // B = A else { - segmentminus=fil->segments[seg-1]; + segmentminus=fil->segments[seg-1+fil->frontseg]; segment->xyzfront[0]=segmentminus->xyzback[0]; // x_i = prior end segment->xyzfront[1]=segmentminus->xyzback[1]; segment->xyzfront[2]=segmentminus->xyzback[2]; Sph_DcmxDcm(segment->dcm,segmentminus->adcm,segment->adcm); } // B_i = A_i . B_{i-1} Sph_DcmtxUnit(segment->adcm,'x',segment->xyzback,segment->xyzfront,segment->len); // x_{i+1} = x_i + l_i * BT_i . xhat - fil->nbs++; } + fil->nseg++; } else { - if(fil->frontbs==0) filArrayShift(fil,0); - if(fil->frontbs==0) filArrayShift(fil,1); // used if nbs=maxbs-1 - seg=fil->frontbs; - segment=fil->segments[seg]; + if(fil->frontseg==0) filArrayShift(fil,0); + if(fil->frontseg==0) filArrayShift(fil,1); // used if nseg=maxseg-1 + seg=0; + segment=fil->segments[seg+fil->frontseg]; segment->len=length; segment->thk=thickness; - if(fil->nbs==0) { + segment->index=seg; + if(fil->nseg==0) { Sph_Xyz2Dcmt(angle,segment->adcm); // B_0 = Dcmt(angle) segment->xyzback[0]=x[0]; // back of segment = input value segment->xyzback[1]=x[1]; segment->xyzback[2]=x[2]; } else { - segmentplus=fil->segments[seg+1]; + segmentplus=fil->segments[seg+1+fil->frontseg]; segment->xyzback[0]=segmentplus->xyzfront[0]; // back of segment = next front segment->xyzback[1]=segmentplus->xyzfront[1]; segment->xyzback[2]=segmentplus->xyzfront[2]; @@ -874,8 +814,9 @@ int filAddSegment(filamentptr fil,const double *x,double length,const double *an Sph_Dcm2Dcm(segment->adcm,segment->dcm); // A_i = B_i Sph_Dcm2Xyz(segment->dcm,segment->ypr); // a_0 = Xyz(B_0) Sph_DcmtxUnit(segment->adcm,'x',segment->xyzfront,segment->xyzback,-segment->len); // x_i = x_{i+1} - l_i * BT_i . xhat - fil->frontbs--; - fil->nbs++; } + fil->frontseg--; + fil->nseg++; + filUpdateSegmentIndices(fil); } return 0; } @@ -896,7 +837,7 @@ int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char varvalues=sim->varvalues; nvar=sim->nvar; - if(fil->nbs==0) { // new filament + if(fil->nseg==0) { // new filament systemrandpos(sim,pos); if(!strcmp(xstr,"u")); else if(strmathsscanf(xstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[0]=f1; @@ -909,16 +850,18 @@ int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char else if(strmathsscanf(zstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[2]=f1; else return 2; - angle[0]=(dim==2 ? unirandCOD(0,2*PI) : thetarandCCD()); + angle[0]=(dim==2 ? unirandCOD(-PI,PI) : thetarandCCD()); angle[1]=unirandCOD(0,2*PI); angle[2]=unirandCOD(0,2*PI); if(!strcmp(thtstr,"u")); else if(strmathsscanf(thtstr,"%mlg",varnames,varvalues,nvar,&f1)==1) angle[0]=f1; else return 3; - if(!strcmp(phistr,"u")); + if(dim==2) angle[1]=0; + else if(!strcmp(phistr,"u")); else if(strmathsscanf(phistr,"%mlg",varnames,varvalues,nvar,&f1)==1) angle[1]=f1; else return 3; - if(!strcmp(chistr,"u")); + if(dim==2) angle[2]=0; + else if(!strcmp(chistr,"u")); else if(strmathsscanf(chistr,"%mlg",varnames,varvalues,nvar,&f1)==1) angle[2]=f1; else return 3; } else { @@ -927,70 +870,31 @@ int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char for(i=0;i1) filRandomAngle(filtype,thickness,angle,1); - if(dim==2) angle[1]=angle[2]=0; + if(fil->nseg>0) { + filRandomAngle(filtype,dim,fil->nseg,thickness,1,angle); + if(dim==2) angle[1]=angle[2]=0; } er=filAddSegment(fil,pos,len,angle,thickness,'b'); if(er) return er; } return 0; } -/* filAddRandomBeads */ -int filAddRandomBeads(filamentptr fil,int number,const char *xstr,const char *ystr,const char *zstr) { - int i,dim,er; - double f1,pos[3],len; - char **varnames; - double *varvalues; - int nvar; - filamenttypeptr filtype; - simptr sim; - - filtype=fil->filtype; - sim=filtype->filss->sim; - dim=sim->dim; - varnames=sim->varnames; - varvalues=sim->varvalues; - nvar=sim->nvar; - - if(fil->nbs==0) { - systemrandpos(sim,pos); - if(!strcmp(xstr,"u")); - else if(strmathsscanf(xstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[0]=f1; - else return 2; - if(!strcmp(ystr,"u")); - else if(strmathsscanf(ystr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[1]=f1; - else return 2; - if(dim==2) pos[2]=0; - else if(!strcmp(zstr,"u")); - else if(strmathsscanf(zstr,"%mlg",varnames,varvalues,nvar,&f1)==1) pos[2]=f1; - else return 2; } - else - pos[0]=pos[1]=pos[2]=0; - - for(i=0;inbs==0) return -1; + if(fil->nseg==0) return -1; if(endchar=='b') - fil->nbs--; + fil->nseg--; else { - seg=++fil->frontbs; - fil->nbs--; - if(!fil->filtype->isbead) { - segment=fil->segments[seg]; - Sph_Dcm2Dcm(segment->adcm,segment->dcm); - Sph_Dcm2Xyz(segment->dcm,segment->ypr); }} + seg=++fil->frontseg; + fil->nseg--; + segment=fil->segments[seg]; + Sph_Dcm2Dcm(segment->adcm,segment->dcm); + Sph_Dcm2Xyz(segment->dcm,segment->ypr); + filUpdateSegmentIndices(fil); } return 0; } @@ -999,20 +903,13 @@ void filTranslate(filamentptr fil,const double *vect,char func) { int seg; double shift[3]; segmentptr segment; - beadptr bead; - seg=fil->frontbs; + seg=fil->frontseg; if(func=='=') { - if(fil->filtype->isbead){ - bead=fil->beads[seg]; - shift[0]=bead->xyz[0]-vect[0]; - shift[1]=bead->xyz[1]-vect[1]; - shift[2]=bead->xyz[2]-vect[2]; } - else{ - segment=fil->segments[seg]; - shift[0]=segment->xyzfront[0]-vect[0]; - shift[1]=segment->xyzfront[1]-vect[1]; - shift[2]=segment->xyzfront[2]-vect[2];}} + segment=fil->segments[seg]; + shift[0]=segment->xyzfront[0]-vect[0]; + shift[1]=segment->xyzfront[1]-vect[1]; + shift[2]=segment->xyzfront[2]-vect[2];} else if(func=='-') { shift[0]=-vect[0]; shift[1]=-vect[1]; @@ -1022,94 +919,142 @@ void filTranslate(filamentptr fil,const double *vect,char func) { shift[1]=vect[1]; shift[2]=vect[2]; } - if(fil->filtype->isbead){ - for(seg=0;segnbs;seg++) { - bead = fil->beads[seg+fil->frontbs]; - bead->xyz[0]+=shift[0]; - bead->xyz[1]+=shift[1]; - bead->xyz[2]+=shift[2]; - bead->xyzold[0]+=shift[0]; - bead->xyzold[1]+=shift[1]; - bead->xyzold[2]+=shift[2]; }} - else { - for(seg=0;segnbs;seg++) { - segment=fil->segments[seg+fil->frontbs]; - segment->xyzfront[0]+=shift[0]; - segment->xyzfront[1]+=shift[1]; - segment->xyzfront[2]+=shift[2]; - segment->xyzback[0]+=shift[0]; - segment->xyzback[1]+=shift[1]; - segment->xyzback[2]+=shift[2]; }} + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg+fil->frontseg]; + segment->xyzfront[0]+=shift[0]; + segment->xyzfront[1]+=shift[1]; + segment->xyzfront[2]+=shift[2]; + segment->xyzback[0]+=shift[0]; + segment->xyzback[1]+=shift[1]; + segment->xyzback[2]+=shift[2]; } return; } /* filRotateVertex */ -void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func) { - (void)fil; - (void)seg; - (void)angle; - (void)endchar; - (void)func; +void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func) { //?? Not written yet + int i; + double dcmdelta[9]; + segmentptr segment,segmentm1,segmentp1; + + segment=fil->segments[fil->frontseg+seg]; + Sph_Xyz2Dcm(angle,dcmdelta); + if(func=='=') Sph_Dcm2Dcm(dcmdelta,segment->dcm); + else if(func=='+') Sph_DcmxDcm(dcmdelta,segment->dcm,segment->dcm); + else Sph_DcmtxDcm(dcmdelta,segment->dcm,segment->dcm); + + if(endchar=='b') { + for(i=seg;inseg;i++) { + segment=fil->segments[fil->frontseg+i]; + if(i==0) Sph_Dcm2Dcm(segment->dcm,segment->adcm); + else { + segmentm1=fil->segments[fil->frontseg+i-1]; + Sph_DcmxDcm(segment->dcm,segmentm1->adcm,segment->adcm); + segment->xyzfront[0]=segmentm1->xyzback[0]; + segment->xyzfront[1]=segmentm1->xyzback[1]; + segment->xyzfront[2]=segmentm1->xyzback[2]; } + Sph_Dcm2Xyz(segment->dcm,segment->ypr); + segment->xyzback[0]=segment->xyzfront[0]+segment->len*segment->adcm[0]; + segment->xyzback[1]=segment->xyzfront[1]+segment->len*segment->adcm[1]; + segment->xyzback[2]=segment->xyzfront[2]+segment->len*segment->adcm[2]; }} + else { + for(i=seg;i>=0;i--) { + segment=fil->segments[fil->frontseg+i]; + if(i==fil->nseg-1); + else { + segmentp1=fil->segments[fil->frontseg+i+1]; + Sph_DcmtxDcm(segmentp1->dcm,segmentp1->adcm,segment->adcm); + segment->xyzback[0]=segmentp1->xyzfront[0]; + segment->xyzback[1]=segmentp1->xyzfront[1]; + segment->xyzback[2]=segmentp1->xyzfront[2]; } + Sph_Dcm2Xyz(segment->dcm,segment->ypr); + segment->xyzfront[0]=segment->xyzback[0]-segment->len*segment->adcm[0]; + segment->xyzfront[1]=segment->xyzback[1]-segment->len*segment->adcm[1]; + segment->xyzfront[2]=segment->xyzback[2]-segment->len*segment->adcm[2]; }} + return; } /* filLengthenSegment */ void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func) { - (void)fil; - (void)seg; - (void)length; - (void)endchar; - (void)func; + int i; + double lenold,lendelta,xdelta[3]; + segmentptr segment,segment2; + + segment=fil->segments[fil->frontseg+seg]; + lenold=segment->len; + if(func=='=') lendelta=length-lenold; + else if(func=='+') lendelta=length; + else lendelta=lenold-length; + + xdelta[0]=lendelta*segment->adcm[0]; // same as SphDcmtxUnit(segment->adcm,'x',xdelta,NULL,lendelta); + xdelta[1]=lendelta*segment->adcm[1]; + xdelta[2]=lendelta*segment->adcm[2]; + if(endchar=='b') { + segment->xyzback[0]+=xdelta[0]; + segment->xyzback[1]+=xdelta[1]; + segment->xyzback[2]+=xdelta[2]; + for(i=seg+1;inseg;i++) { + segment2=fil->segments[fil->frontseg+i]; + segment2->xyzfront[0]=segment->xyzback[0]; + segment2->xyzfront[1]=segment->xyzback[1]; + segment2->xyzfront[2]=segment->xyzback[2]; + segment2->xyzback[0]+=xdelta[0]; + segment2->xyzback[1]+=xdelta[1]; + segment2->xyzback[2]+=xdelta[2]; + segment=segment2; }} + else { + segment->xyzfront[0]-=xdelta[0]; + segment->xyzfront[1]-=xdelta[1]; + segment->xyzfront[2]-=xdelta[2]; + for(i=seg-1;i>=0;i--) { + segment2=fil->segments[fil->frontseg+i]; + segment2->xyzback[0]=segment->xyzfront[0]; + segment2->xyzback[1]=segment->xyzfront[1]; + segment2->xyzback[2]=segment->xyzfront[2]; + segment2->xyzfront[0]-=xdelta[0]; + segment2->xyzfront[1]-=xdelta[1]; + segment2->xyzfront[2]-=xdelta[2]; + segment=segment2; }} + return; } /* filReverseFilament */ -void filReverseFilament(filamentptr fil) { - (void)fil; +void filReverseFilament(filamentptr fil) { //?? Not written yet + filUpdateSegmentIndices(fil); return; } /* filCopyFilament */ int filCopyFilament(filamentptr filto,const filamentptr filfrom) { - int isbead,i; - beadptr beadto,beadfrom; + int i; segmentptr segmentto,segmentfrom; if(!filto || !filfrom) return 2; - if(filfrom->filtype) isbead=filfrom->filtype->isbead; - else if(filfrom->beads) isbead=1; - else isbead=0; - filto->filtype=filfrom->filtype; - filto->nbs=0; - filto->frontbs=0; + filto->nseg=0; + filto->frontseg=0; filto->nbranch=0; filto->nmonomer=0; filto->frontmonomer=0; - filto=filalloc(filto,isbead,filfrom->nbs,filfrom->nbranch,filfrom->nmonomer); + filto=filalloc(filto,filfrom->nseg,filfrom->nbranch,filfrom->nmonomer); if(!filto) return 1; - if(isbead) { - for(i=0;inbs;i++) { - beadto=filto->beads[i]; - beadfrom=filfrom->beads[i+filfrom->frontbs]; - copyVD(beadfrom->xyz,beadto->xyz,3); - copyVD(beadfrom->xyzold,beadto->xyzold,3); }} - else { - for(i=0;inbs;i++) { - segmentto=filto->segments[i]; - segmentfrom=filfrom->segments[i+filfrom->frontbs]; - copyVD(segmentfrom->xyzfront,segmentto->xyzfront,3); - copyVD(segmentfrom->xyzback,segmentto->xyzback,3); - segmentto->len=segmentfrom->len; - copyVD(segmentfrom->ypr,segmentto->ypr,3); - copyVD(segmentfrom->dcm,segmentto->dcm,9); - copyVD(segmentfrom->adcm,segmentto->adcm,9); - segmentto->thk=segmentfrom->thk; }} - filto->nbs=filfrom->nbs; + for(i=0;inseg;i++) { + segmentto=filto->segments[i]; + segmentfrom=filfrom->segments[i+filfrom->frontseg]; + segmentto->index=segmentfrom->index; + copyVD(segmentfrom->xyzfront,segmentto->xyzfront,3); + copyVD(segmentfrom->xyzback,segmentto->xyzback,3); + segmentto->len=segmentfrom->len; + segmentto->thk=segmentfrom->thk; + copyVD(segmentfrom->ypr,segmentto->ypr,3); + copyVD(segmentfrom->dcm,segmentto->dcm,9); + copyVD(segmentfrom->adcm,segmentto->adcm,9); } + filto->nseg=filfrom->nseg; filto->frontend=filfrom->frontend; filto->backend=filfrom->backend; @@ -1125,6 +1070,31 @@ int filCopyFilament(filamentptr filto,const filamentptr filfrom) { return 0; } +/* filAddFilament */ +filamentptr filAddFilament(filamenttypeptr filtype,const char *filname) { + int f; + filamentptr fil; + char autoname[STRCHAR]; + + if(filname) { + f=stringfind(filtype->filnames,filtype->nfil,filname); + if(f>=0) return filtype->fillist[f]; } + else { + sprintf(autoname,"%s.%i",filtype->ftname,filtype->autonamenum++); + filname=autoname; } + + if(filtype->nfil==filtype->maxfil) { + filtype=filamenttypealloc(filtype,filtype->maxfil*2+1,0); + if(!filtype) return NULL; } + f=filtype->nfil++; + strncpy(filtype->filnames[f],filname,STRCHAR-1); + filtype->filnames[f][STRCHAR-1]='\0'; + fil=filtype->fillist[f]; + filsetcondition(filtype->filss,SClists,0); + + return fil; } + + /******************************************************************************/ /********************************* filament type ******************************/ /******************************************************************************/ @@ -1136,36 +1106,45 @@ int filtypeSetParam(filamenttypeptr filtype,const char *param,int index,double v er=0; if(!strcmp(param,"stdlen")) { - if(value<=0) er=2; + if(value<0) er=2; else filtype->stdlen=value; } else if(!strcmp(param,"stdypr")) { - if(index<0) + if(value<-PI || value>PI) er=2; + else if(index<0) filtype->stdypr[0]=filtype->stdypr[1]=filtype->stdypr[2]=value; else filtype->stdypr[index]=value; } - else if(!strcmp(param,"klen")) - filtype->klen=value; + else if(!strcmp(param,"klen")) { + filtype->klen=value; } else if(!strcmp(param,"kypr")) { if(index<0) filtype->kypr[0]=filtype->kypr[1]=filtype->kypr[2]=value; else filtype->kypr[index]=value; } else if(!strcmp(param,"kT")) { - if(value<=0) er=2; + if(value<0) er=2; else filtype->kT=value; } - else if(!strcmp(param,"treadrate")) - filtype->treadrate=value; + else if(!strcmp(param,"treadrate")) { + if(value<0) er=2; + else filtype->treadrate=value; } else if(!strcmp(param,"viscosity")) { if(value<=0) er=2; else filtype->viscosity=value; } - else if(!strcmp(param,"beadradius")) { + else if(!strcmp(param,"bundle")) { + if(value<=0) er=2; + else filtype->bundlevalue=value; } + + else if(!strcmp(param,"radius")) { if(value<=0) er=2; else filtype->filradius=value; } + else if(!strcmp(param,"facetwist")) { + filtype->facetwist=value; } + return er; } @@ -1215,12 +1194,7 @@ int filtypeSetShiny(filamenttypeptr filtype,double shiny) { /* filtypeSetDynamics */ int filtypeSetDynamics(filamenttypeptr filtype,enum FilamentDynamics fd) { - filtype->dynamics=fd; //?? check for pre-existing dynamics and convert if needed - if(fd==FDrigidbeads) filtype->isbead=1; - else if(fd==FDrigidsegments) filtype->isbead=0; - else if(fd==FDrouse) filtype->isbead=1; - else if(fd==FDalberts) filtype->isbead=0; - else if(fd==FDnedelec) filtype->isbead=0; + filtype->dynamics=fd; return 0; } @@ -1230,6 +1204,43 @@ int filtypeSetBiology(filamenttypeptr filtype,enum FilamentBiology fb) { return 0; } +/* filtypeAddFace */ +int filtypeAddFace(filamenttypeptr filtype,const char* facename) { + if(filtype->nface==filtype->maxface) { + filtype=filamenttypealloc(filtype,0,filtype->maxface*2+1); + if(!filtype) return -1; } + strcpy(filtype->facename[filtype->nface++],facename); + return 0; } + + +/* filAddFilamentType */ +filamenttypeptr filAddFilamentType(simptr sim,const char *ftname) { + int ft,er; + filamentssptr filss; + filamenttypeptr filtype; + + if(!sim->filss) { + er=filenablefilaments(sim); + if(er) return NULL; } + filss=sim->filss; + + ft=stringfind(filss->ftnames,filss->ntype,ftname); + if(ft<0) { + if(filss->ntype==filss->maxtype) { + filss=filssalloc(filss,filss->maxtype*2+1); + if(!filss) return NULL; } + ft=filss->ntype++; + strncpy(filss->ftnames[ft],ftname,STRCHAR-1); + filss->ftnames[ft][STRCHAR-1]='\0'; + filtype=filss->filtypes[ft]; + filsetcondition(filss,SClists,0); } + else + filtype=filss->filtypes[ft]; + + return filtype; } + + + /******************************************************************************/ /**************************** filament superstructure *************************/ /******************************************************************************/ @@ -1253,7 +1264,7 @@ int filenablefilaments(simptr sim) { if(sim->filss) return 0; - filss=filssalloc(sim->filss,1); + filss=filssalloc(sim->filss,0); if(!filss) return 1; sim->filss=filss; filss->sim=sim; @@ -1261,66 +1272,34 @@ int filenablefilaments(simptr sim) { return 0; } -/* filAddFilament */ -filamentptr filAddFilament(filamenttypeptr filtype,filamentptr fil,const char *filname) { - int f; - filamentptr fil2; - - if(!filtype && fil) // no filtype, yes fil - return fil; - - else if(!filtype) { // no filtype, no fil - fil=filalloc(NULL,0,0,0,0); - if(!fil) return NULL; - fil->filname=EmptyString(); - if(!fil->filname) return NULL; } +/* filupdateparams */ +int filupdateparams(simptr sim) { + (void) sim; + return 0; } - else { // yes filtype, yes or no fil - f=stringfind(filtype->filnames,filtype->nfil,filname); - if(f<0) { - if(filtype->nfil==filtype->maxfil) { - filtype=filamenttypealloc(filtype,filtype->maxfil*2+1,0); - if(!filtype) return NULL; } - f=filtype->nfil++; - strncpy(filtype->filnames[f],filname,STRCHAR-1); - filtype->filnames[f][STRCHAR-1]='\0'; - fil2=filtype->fillist[f]; - if(fil) { - filCopyFilament(fil2,fil); - fil2->filtype=filtype; - free(fil->filname); - filfree(fil); } - fil=fil2; - filsetcondition(filtype->filss,SClists,0); } - else - fil=filtype->fillist[f]; } - return fil; } +/* filupdatelists */ +int filupdatelists(simptr sim) { + (void) sim; + return 0; } -/* filAddFilamentType */ -filamenttypeptr filAddFilamentType(simptr sim,const char *ftname) { - int ft; +/* filupdate */ +int filupdate(simptr sim) { + int er; filamentssptr filss; - filamenttypeptr filtype; - filamentssptr fssptr; filss=sim->filss; - - ft=stringfind(filss->ftnames,filss->ntype,ftname); - if(ft<0) { - if(filss->ntype==filss->maxtype) { - fssptr=filssalloc(filss,filss->maxtype*2+1); - if(!fssptr) return NULL; } - ft=filss->ntype++; - strncpy(filss->ftnames[ft],ftname,STRCHAR-1); - filss->ftnames[ft][STRCHAR-1]='\0'; - filtype=filss->filtypes[ft]; - filsetcondition(filss,SClists,0); } - else - filtype=filss->filtypes[ft]; - - return filtype; } + if(filss) { + if(filss->condition<=SClists) { + er=filupdatelists(sim); + if(er) return er; + filsetcondition(filss,SCparams,1); } + if(filss->condition==SCparams) { + er=filupdateparams(sim); + if(er) return er; + filsetcondition(filss,SCok,1); }} + return 0; } /******************************************************************************/ @@ -1332,10 +1311,8 @@ filamenttypeptr filAddFilamentType(simptr sim,const char *ftname) { filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr filtype,const char *word,char *line2) { char **varnames; double *varvalues; - int nvar; + int nvar,dim; - filamentssptr filss; - UNUSED(filss); int itct,er,i1,i2; char nm[STRCHARLONG],nm1[STRCHARLONG]; double fltv1[9],f1; @@ -1344,7 +1321,7 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi enum FilamentBiology fb; printf("%s\n",word);//?? debug - filss=sim->filss; + dim=sim->dim; varnames=sim->varnames; varvalues=sim->varvalues; @@ -1360,10 +1337,10 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi else if(!strcmp(word,"dynamics")) { // dynamics CHECKS(filtype,"need to enter filament type name before dynamics"); itct=sscanf(line2,"%s",nm1); - CHECKS(itct==1,"dynamics options: RigidBeads, RigidSegments, Rouse, Alberts, Nedelec"); + CHECKS(itct==1,"dynamics options: none, newton"); fd=filstring2FD(nm1); - CHECKS(fd!=FDnone,"dynamics options: RigidBeads, RigidSegments, Rouse, Alberts, Nedelec"); - er=filtypeSetDynamics(filtype,fd); //?? beware of changing between beads and segments + CHECKS(fd!=FDnone,"dynamics options: none, newton"); + er=filtypeSetDynamics(filtype,fd); CHECKS(!er,"BUG: error setting filament dynamics"); CHECKS(!strnword(line2,2),"unexpected text following dynamics"); } @@ -1429,7 +1406,7 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi CHECKS(filtype,"need to enter filament type name before kT"); itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1); CHECKS(itct==1,"kT format: value"); - CHECKS(f1>0,"kT value needs to be greater than 0"); + CHECKS(f1>=0,"kT value needs to be >=0"); filtypeSetParam(filtype,"kT",0,f1); CHECKS(!strnword(line2,2),"unexpected text following kT"); } @@ -1437,8 +1414,9 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi CHECKS(filtype,"need to enter filament type name before treadmill_rate"); itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1); CHECKS(itct==1,"treadmill_rate format: value"); + CHECKS(f1>=0,"treadmill_rate cannot be negative"); filtypeSetParam(filtype,"treadrate",0,f1); - CHECKS(!strnword(line2,2),"unexpected text following treadrate"); } + CHECKS(!strnword(line2,2),"unexpected text following treadmill_rate"); } else if(!strcmp(word,"viscosity")) { // viscosity CHECKS(filtype,"need to enter filament type name before viscosity"); @@ -1448,29 +1426,27 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi filtypeSetParam(filtype,"viscosity",0,f1); CHECKS(!strnword(line2,2),"unexpected text following viscosity"); } - else if(!strcmp(word,"bead_radius")) { // beadradius - CHECKS(filtype,"need to enter filament type name before bead radius"); - itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1); - CHECKS(itct==1,"bead_radius format: value"); - filtypeSetParam(filtype,"beadradius",0,f1); - CHECKS(!strnword(line2,2),"unexpected text following bead radius"); } - else if(!strcmp(word,"standard_length")) { // standard_length CHECKS(filtype,"need to enter filament type name before standard_length"); itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&f1); CHECKS(itct==1,"standard_length format: value"); - CHECKS(f1>0,"standard_length value needs to be greater than 0"); + CHECKS(f1>=0,"standard_length value needs to be >=0"); filtypeSetParam(filtype,"stdlen",0,f1); CHECKS(!strnword(line2,2),"unexpected text following standard_length"); } else if(!strcmp(word,"standard_angle")) { // standard_angle CHECKS(filtype,"need to enter filament type name before standard_angle"); - itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]); - CHECKS(itct==3,"standard_angle format: yaw pitch roll"); + if(dim==2) { + itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&fltv1[0]); + CHECKS(itct==1,"standard_angle format: angle"); + fltv1[1]=fltv1[2]=0; } + else { + itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]); + CHECKS(itct==3,"standard_angle format: yaw pitch roll"); } filtypeSetParam(filtype,"stdypr",0,fltv1[0]); filtypeSetParam(filtype,"stdypr",1,fltv1[1]); filtypeSetParam(filtype,"stdypr",2,fltv1[2]); - CHECKS(!strnword(line2,4),"unexpected text following standard_angle"); } + CHECKS(!strnword(line2,dim==2?2:4),"unexpected text following standard_angle"); } else if(!strcmp(word,"force_length")) { // force_length CHECKS(filtype,"need to enter filament type name before force_length"); @@ -1481,8 +1457,13 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi else if(!strcmp(word,"force_angle")) { // force_angle CHECKS(filtype,"need to enter filament type name before force_angle"); - itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]); - CHECKS(itct==3,"force_angle format: yaw pitch roll"); + if(dim==2) { + itct=strmathsscanf(line2,"%mlg",varnames,varvalues,nvar,&fltv1[0]); + CHECKS(itct==1,"force_angle format: value"); + fltv1[1]=fltv1[2]=-1; } + else { + itct=strmathsscanf(line2,"%mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2]); + CHECKS(itct==3,"force_angle format: yaw pitch roll"); } filtypeSetParam(filtype,"kypr",0,fltv1[0]); filtypeSetParam(filtype,"kypr",1,fltv1[1]); filtypeSetParam(filtype,"kypr",2,fltv1[2]); @@ -1519,28 +1500,19 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty varvalues=sim->varvalues; nvar=sim->nvar; - if(!strcmp(word,"name")) { // name - itct=sscanf(line2,"%s",nm); - CHECKS(itct==1,"error reading filament name"); - fil=filAddFilament(filtype,fil,nm); - CHECKS(fil,"failed to add filament"); - CHECKS(!strnword(line2,2),"unexpected text following filament name"); } - - else if(!strcmp(word,"type")) { // type - itct=sscanf(line2,"%s",nm1); - CHECKS(itct==1,"error reading filament type"); + if(!strcmp(word,"name")) { + itct=sscanf(line2,"%s %s",nm1,nm); + CHECKS(itct==2,"correct format: filament_type filament_name"); i1=stringfind(filss->ftnames,filss->ntype,nm1); CHECKS(i1>=0,"filament type not recognized"); filtype=filss->filtypes[i1]; - if(fil) { - CHECKS(fil->filtype==NULL || fil->filtype==filtype,"filament type was already defined and cannot be changed"); - fil=filAddFilament(filtype,fil,fil->filname); - CHECKS(fil,"failed to add filament"); } - CHECKS(!strnword(line2,2),"unexpected text following name"); } + fil=filAddFilament(filtype,nm); + CHECKS(fil,"failed to add filament"); + CHECKS(!strnword(line2,3),"unexpected text following filament name"); } else if(!strcmp(word,"first_segment")) { // first_segment CHECKS(fil,"need to enter filament name before first_segment"); - CHECKS(fil->nbs==0,"filament already has segments in it"); + CHECKS(fil->nseg==0,"filament already has segments in it"); itct=strmathsscanf(line2,"%mlg %mlg %mlg %mlg %mlg %mlg %mlg",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2],&length,&angle[0],&angle[1],&angle[2]); CHECKS(itct==7,"first_segment format: x y z length angle0 angle1 angle2 [thickness]"); CHECKS(length>0,"length needs to be >0"); @@ -1557,7 +1529,7 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty else if(!strcmp(word,"add_segment")) { // add_segment CHECKS(fil,"need to enter filament name before add_segment"); - CHECKS(fil->nbs>0,"use first_segment to enter the first segment"); + CHECKS(fil->nseg>0,"use first_segment to enter the first segment"); itct=strmathsscanf(line2,"%mlg %mlg %mlg %mlg",varnames,varvalues,nvar,&length,&angle[0],&angle[1],&angle[2]); CHECKS(itct==4,"add_segment format: length angle0 angle1 angle2 [thickness [end]]"); CHECKS(length>0,"length needs to be >0"); @@ -1582,7 +1554,7 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty else if(!strcmp(word,"remove_segment")) { // remove_segment CHECKS(fil,"need to enter filament name before remove_segment"); - CHECKS(fil->nbs>0,"filament has no segments to remove"); + CHECKS(fil->nseg>0,"filament has no segments to remove"); itct=sscanf(line2,"%s",nm1); CHECKS(itct==1,"remove_segment format: end"); endchar='b'; @@ -1596,11 +1568,12 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty else if(!strcmp(word,"random_segments")) { // random_segments CHECKS(fil,"need to enter filament name before random_segments"); CHECKS(fil->filtype,"need to enter filament type before random_segments"); + CHECKS(fil->filtype->klen==-1 || fil->filtype->klen>0,"cannot compute random segments because the filament type has length force constant equals 0"); itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1); CHECKS(itct==1,"random_segments format: number [x y z theta phi chi] [thickness]"); CHECKS(i1>0,"number needs to be >0"); line2=strnword(line2,2); - if(fil->nbs==0) { + if(fil->nseg==0) { CHECKS(line2,"missing position and angle information"); itct=sscanf(line2,"%s %s %s %s %s %s",str1,str2,str3,str4,str5,str6); CHECKS(itct==6,"random_segments format: number [x y z theta phi chi] [thickness]"); @@ -1636,7 +1609,7 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty CHECKS(fil,"need to enter filament name before copy"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"error reading filament name"); - fil2=filAddFilament(fil->filtype,NULL,nm); + fil2=filAddFilament(fil->filtype,nm); CHECKS(fil2,"failed to add filament"); er=filCopyFilament(fil2,fil); CHECKS(!strnword(line2,2),"unexpected text following copy"); } @@ -1698,17 +1671,19 @@ int filloadtype(simptr sim,ParseFilePtr *pfpptr,char *line2) { /* filloadfil */ -int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2,filamenttypeptr filtype) { +int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2) { ParseFilePtr pfp; char word[STRCHAR],errstring[STRCHARLONG]; int done,pfpcode,firstline2; filamentptr fil; + filamenttypeptr filtype; pfp=*pfpptr; done=0; firstline2=line2?1:0; filenablefilaments(sim); fil=NULL; + filtype=NULL; while(!done) { if(pfp->lctr==0) @@ -1720,21 +1695,21 @@ int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2,filamenttypeptr filty else pfpcode=Parse_ReadLine(&pfp,word,&line2,errstring); - *pfpptr=pfp; - CHECKS(pfpcode!=3,"%s",errstring); - - if(pfpcode==0); // already taken care of - else if(pfpcode==2) { // end reading - done=1; } - else if(!strcmp(word,"end_filament")) { // end_filament_type - CHECKS(!line2,"unexpected text following end_filament"); - return 0; } - else if(!line2) { // just word - CHECKS(0,"unknown word or missing parameter"); } - else { - fil=filreadstring(sim,pfp,fil,filtype,word,line2); - CHECK(fil); // failed but error has already been sent - filtype=fil->filtype; }} + *pfpptr=pfp; + CHECKS(pfpcode!=3,"%s",errstring); + + if(pfpcode==0); // already taken care of + else if(pfpcode==2) { // end reading + done=1; } + else if(!strcmp(word,"end_filament")) { // end_filament_type + CHECKS(!line2,"unexpected text following end_filament"); + return 0; } + else if(!line2) { // just word + CHECKS(0,"unknown word or missing parameter"); } + else { + fil=filreadstring(sim,pfp,fil,filtype,word,line2); + CHECK(fil); // failed but error has already been sent + filtype=fil->filtype; }} CHECKS(0,"end of file encountered before end_filament statement"); // end of file @@ -1744,167 +1719,112 @@ int filloadfil(simptr sim,ParseFilePtr *pfpptr,char *line2,filamenttypeptr filty return 1; } -/******************************************************************************/ -/******************************* structure updating ***************************/ -/******************************************************************************/ - - -/* filupdateparams */ -int filupdateparams(simptr sim) { - (void) sim; - return 0; } - - -/* filupdatelists */ -int filupdatelists(simptr sim) { - (void) sim; - return 0; } - - -/* filupdate */ -int filsupdate(simptr sim) { - int er; - filamentssptr filss; - - filss=sim->filss; - if(filss) { - if(filss->condition<=SClists) { - er=filupdatelists(sim); - if(er) return er; - filsetcondition(filss,SCparams,1); } - if(filss->condition==SCparams) { - er=filupdateparams(sim); - if(er) return er; - filsetcondition(filss,SCok,1); }} - return 0; } - - /******************************************************************************/ /**************************** core simulation functions ***********************/ /******************************************************************************/ /* filSegmentXSurface */ -int filSegmentXSurface(simptr sim,filamentptr fil,char endchar) { //?? Currently only suitable for treadmilling, not dynamics - int s,p,mxs; +int filSegmentXSurface(simptr sim,segmentptr segment,panelptr *pnlptr) { + int s,p,cross; surfaceptr srf; panelptr pnl; double *pt1,*pt2,crosspt[3]; enum PanelShape ps; - filamenttypeptr filtype; - beadptr bead,beadplus; - segmentptr segment; - - filtype = fil->filtype; if(!sim->srfss) return 0; - if(endchar=='f') { - if(filtype->isbead){ - bead=fil->beads[fil->frontbs]; - beadplus=fil->beads[fil->frontbs+1]; - pt1=bead->xyz; - pt2=beadplus->xyz; } - else{ - segment=fil->segments[fil->frontbs]; - pt1=segment->xyzfront; - pt2=segment->xyzback; }} - else { - if(filtype->isbead){ - bead=fil->beads[fil->nbs+fil->frontbs-1]; - beadplus=fil->beads[fil->nbs+fil->frontbs]; - pt1=bead->xyz; - pt2=beadplus->xyz; } - else{ - segment=fil->segments[fil->nbs+fil->frontbs]; - pt1=segment->xyzfront; - pt2=segment->xyzback; }} - - mxs=0; - for(s=0;ssrfss->nsrf && !mxs;s++) { + pt1=segment->xyzfront; + pt2=segment->xyzback; + + cross=0; + for(s=0;ssrfss->nsrf && !cross;s++) { srf=sim->srfss->srflist[s]; - for(ps=(enum PanelShape)0;psnpanel[ps] && !mxs;p++) { + for(ps=(enum PanelShape)0;psnpanel[ps] && !cross;p++) { pnl=srf->panels[ps][p]; - mxs=lineXpanel(pt1,pt2,pnl,3,crosspt,NULL,NULL,NULL,NULL,NULL,0); }} -// printf("pt1=(%g %g %g), pt2=(%g %g %g), mxs=%i\n",pt1[0],pt1[1],pt1[2],pt2[0],pt2[1],pt2[2],mxs); }} - - return mxs; } + cross=lineXpanel(pt1,pt2,pnl,sim->dim,crosspt,NULL,NULL,NULL,NULL,NULL,0); }} + if(cross && pnlptr) *pnlptr=pnl; + return cross; } /* filSegmentXFilament */ -int filSegmentXFilament(simptr sim,filamentptr fil,char endchar,filamentptr *filptr) { //?? Currently only supported for segments not beads - int f,i,mxf,mn,mn1,ft; - double thk=0.0,*pt1,*pt2,dist=0.0; +int filSegmentXFilament(simptr sim,segmentptr segment,filamentptr *filptr) { + int f,i,ft,cross; + double thk,*ptf,*ptb,dist; filamentssptr filss; filamenttypeptr filtype; - filamentptr fil2; - segmentptr segment; + filamentptr fil,fil2; + segmentptr segmentm1,segmentp1,segment2; - if(endchar=='f') { - segment=fil->segments[fil->frontbs]; - pt1=segment->xyzfront; - pt2=segment->xyzback; - thk=segment->thk; - mn=mn1=fil->frontbs; - if(fil->nbs>1) mn1=fil->frontbs+1; } - else { - segment=fil->segments[fil->nbs+fil->frontbs]; - pt1=segment->xyzfront; - pt2=segment->xyzback; - mn=mn1=fil->nbs+fil->frontbs-1; - if(fil->nbs>1) mn1=fil->nbs+fil->frontbs-2; } + fil=segment->fil; + ptf=segment->xyzfront; + ptb=segment->xyzback; + thk=segment->thk; + segmentm1=(segment->index==0) ? NULL:fil->segments[fil->frontseg+segment->index-1]; + segmentp1=(segment->index==fil->nseg-1) ? NULL:fil->segments[fil->frontseg+segment->index+1]; - mxf=0; + cross=0; filss=sim->filss; - fil2=NULL; - for(ft=0;ftntype && !mxf;ft++) { + for(ft=0;ftntype && !cross;ft++) { filtype=filss->filtypes[ft]; - for(f=0;fnfil && !mxf;f++) { + for(f=0;fnfil && !cross;f++) { fil2=filtype->fillist[f]; - for(i=fil2->frontbs;inbs+fil2->frontbs && !mxf;i++) { - if(!(fil2==fil && (i==mn || i==mn1))) { - dist=Geo_NearestSeg2SegDist(pt1,pt2,fil2->segments[i]->xyzfront,fil2->segments[i]->xyzfront); //?? check - if(distsegments[i]->thk) mxf=1; }}}} - if(mxf && filptr) + for(i=0;inseg && !cross;i++) { + segment2=fil2->segments[i+fil2->frontseg]; + if(!(segment2==segment || segment2==segmentm1 || segment2==segmentp1)) { + dist=Geo_NearestSeg2SegDist(ptf,ptb,segment2->xyzfront,segment2->xyzback); + if(distthk) cross=1; }}}} + if(cross && filptr) *filptr=fil2; - return mxf; } + return cross; } -/* filTreadmill */ -void filTreadmill(simptr sim,filamentptr fil,int steps) { //?? Currently for segments only - int i,mxs,done,iter; - double thk,length,angle[3],sigmamult,angletry[3],dcm[9]; - filamentptr fil2; +/* filAddOneRandomSegment */ +int filAddOneRandomSegment(simptr sim,filamentptr fil,const double *x,double thickness,char endchar,int constraints) { + int dim,er,iter,cross,tryagain; filamenttypeptr filtype; + double len,angle[3]; + panelptr pnl; + segmentptr segment; filtype=fil->filtype; - for(i=0;isegments[0]->thk; - done=0; - sigmamult=1; - angletry[0]=angletry[1]=angletry[2]=0; - for(iter=0;iter<500 && !done;iter++) { - length=filRandomLength(filtype,thk,sigmamult); - filRandomAngle(filtype,thk,angle,sigmamult); - if(iter>0 && coinrandD(0.5)) - filAddSegment(fil,NULL,length,angletry,thk,'b'); - else - filAddSegment(fil,NULL,length,angle,thk,'b'); - mxs=filSegmentXSurface(sim,fil,'b'); - if(!mxs) { - mxs=filSegmentXFilament(sim,fil,'b',&fil2); - if(mxs) { - mxs=coinrandD(0.995); - Sph_DcmxDcmt(fil2->segments[fil2->nbs+fil2->frontbs-1]->adcm,fil->segments[fil->nbs+fil->frontbs-2]->adcm,dcm); - Sph_Dcm2Xyz(dcm,angletry); }} - if(mxs) { - filRemoveSegment(fil,'b'); - sigmamult*=1.01; } - else done=1; } - if(done) - filRemoveSegment(fil,'f'); } + dim=sim->dim; - return; } + tryagain=1; + len=filRandomLength(filtype,thickness,1); + filRandomAngle(filtype,dim,fil->nseg,thickness,1,angle); + er=filAddSegment(fil,x,len,angle,thickness,endchar); + if(er) return er; + segment=fil->segments[fil->frontseg+((endchar=='f')?0:fil->nseg-1)]; + if(constraints==0) tryagain=0; + for(iter=0;iternseg,thickness,1,angle); + filLengthenSegment(fil,segment->index,len,endchar,'='); + filRotateVertex(fil,segment->index,angle,endchar,'='); + tryagain=1; }}} + + if(tryagain) { + filRemoveSegment(fil,endchar); + return 2; } + + return 0; } + + +/* filTreadmill */ +int filTreadmill(simptr sim,filamentptr fil,char endchar) { + int er; + + er=0; + if(fil->nseg<1) return 2; + er=filAddOneRandomSegment(sim,fil,NULL,fil->segments[fil->frontseg+((endchar=='b')?fil->nseg-1:0)]->thk,endchar,1); + if(!er) + filRemoveSegment(fil,(endchar=='b')?'f':'b'); + return er; } /* filDynamics */ @@ -1912,51 +1832,50 @@ int filDynamics(simptr sim) { filamentssptr filss; filamentptr fil; filamenttypeptr filtype; - beadptr bead,beadplus,beadminus; - int f,b,d,ft; - double k1,k2; // Force & gaussian constants - int dim; + int f,ft,i,treadnum; filss=sim->filss; if(!filss) return 0; - dim=sim->dim; - for(ft=0;ftntype;ft++) { filtype=filss->filtypes[ft]; - for(f=0;fnfil;f++) { - fil=filtype->fillist[f]; - if(filtype->treadrate>0) - filTreadmill(sim,fil,poisrandD(filtype->treadrate*sim->dt)); - + if(filtype->treadrate>0) + for(f=0;fnfil;f++) { + fil=filtype->fillist[f]; + treadnum=poisrandD(filtype->treadrate*sim->dt); + for(i=0;idynamics==FDrouse) { k1 = 3*filtype->kT*sim->dt/(6*PI*filtype->viscosity*filtype->filradius*filtype->stdlen*filtype->stdlen); //?? Double check this is kuhn length squared k2 = sqrt(2*filtype->kT/(6*PI*filtype->viscosity*filtype->filradius)); - for(b=fil->frontbs;b<=fil->nbs+fil->frontbs;b++){ + for(b=fil->frontseg;b<=fil->nseg+fil->frontseg;b++){ for(d=0;dbeads[b]->xyzold[d]=fil->beads[b]->xyz[d];}} //?? PERHAPS include a check for number of segments to be >= 2 - b=fil->frontbs; + b=fil->frontseg; bead=fil->beads[b]; beadplus=fil->beads[b+1]; for(d=0;dxyz[d]=bead->xyzold[d]-k1*(bead->xyzold[d]-beadplus->xyzold[d])+k2*gaussrandD();} - for(b=fil->frontbs+1;bnbs+fil->frontbs;b++){ + for(b=fil->frontseg+1;bnseg+fil->frontseg;b++){ beadminus=fil->beads[b-1]; bead=fil->beads[b]; beadplus=fil->beads[b+1]; for(d=0;dxyz[d]=bead->xyzold[d]-k1*(2*bead->xyzold[d]-beadminus->xyzold[d]-beadplus->xyzold[d])+k2*gaussrandD();}} - b=fil->nbs+fil->frontbs; + b=fil->nseg+fil->frontseg; beadminus=fil->beads[b-1]; bead=fil->beads[b]; for(d=0;dxyz[d]=bead->xyzold[d]-k1*(bead->xyzold[d]-beadminus->xyzold[d])+k2*gaussrandD();} - }}} + }*/ + }} return 0; } diff --git a/source/Smoldyn/smolgraphics.c b/source/Smoldyn/smolgraphics.c index 6bfaf6ac..faaab9f5 100644 --- a/source/Smoldyn/smolgraphics.c +++ b/source/Smoldyn/smolgraphics.c @@ -1123,15 +1123,11 @@ void RenderFilaments(simptr sim) { glPointSize((GLfloat)filtype->edgepts); glBegin(GL_POINTS); } - for(vtx=fil->frontbs;vtxnbs+fil->frontbs;vtx++) { - if(filtype->isbead) - point=fil->beads[vtx]->xyz; - else - point=fil->segments[vtx]->xyzfront; + for(vtx=fil->frontseg;vtxnseg+fil->frontseg;vtx++) { + point=fil->segments[vtx]->xyzfront; glVertex3d((GLdouble)(point[0]),(GLdouble)(point[1]),(GLdouble)(point[2])); } - if(!filtype->isbead) { - point=fil->segments[vtx-1]->xyzback; - glVertex3d((GLdouble)(point[0]),(GLdouble)(point[1]),(GLdouble)(point[2])); } + point=fil->segments[vtx-1]->xyzback; + glVertex3d((GLdouble)(point[0]),(GLdouble)(point[1]),(GLdouble)(point[2])); glEnd(); } else if(drawmode&DMface) { @@ -1141,7 +1137,7 @@ void RenderFilaments(simptr sim) { //glMaterialfv(GL_FRONT,GL_SPECULAR,gl2Double2GLfloat(srf->fcolor,glfvect,4)); //glMaterialfv(GL_BACK,GL_SPECULAR,gl2Double2GLfloat(srf->bcolor,glfvect,4)); glMateriali(GL_FRONT,GL_SHININESS,(GLint)filtype->shiny); } - for(vtx=fil->frontbs;vtxnbs+fil->frontbs;vtx++) + for(vtx=fil->frontseg;vtxnseg+fil->frontseg;vtx++) /*gl2drawtwistprism(fil->px[vtx],fil->px[vtx+1],fil->nface,fil->po[vtx],twist,fil->radius,fil->facecolor)*/; } if(glIsEnabled(GL_LINE_STIPPLE)) glDisable(GL_LINE_STIPPLE); }} diff --git a/source/Smoldyn/smolsim.cpp b/source/Smoldyn/smolsim.cpp index e65f434a..49da7375 100644 --- a/source/Smoldyn/smolsim.cpp +++ b/source/Smoldyn/smolsim.cpp @@ -591,7 +591,7 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { char nm[STRCHAR],nm1[STRCHAR],shapenm[STRCHAR],ch,rname[STRCHAR],fname[STRCHAR],pattern[STRCHAR]; char str1[STRCHAR],str2[STRCHAR],str3[STRCHAR],str4[STRCHAR],str5[STRCHAR],str6[STRCHAR]; char errstr[STRCHARLONG]; - int er,i,nmol,d,i1,s,c,ll,order,*index,ft; + int er,i,nmol,d,i1,s,c,ll,order,*index,ft,f; int rulelist[MAXORDER+MAXPRODUCT],r,ord,rct,prd,itct,prt,lt,detailsi[8]; long int pserno,sernolist[MAXPRODUCT]; double flt1,flt2,v1[DIMMAX*DIMMAX],v2[4],poslo[DIMMAX],poshi[DIMMAX],thick; @@ -1483,9 +1483,8 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { else if(!strcmp(word,"max_filament")) { // max_filament CHECKS(0,"max_filament has been deprecated"); } -/* else if(!strcmp(word,"new_filament")) { // new_filament - fil=filreadstring(sim,pfp,NULL,"name",line2); + fil=filreadstring(sim,pfp,NULL,NULL,"name",line2); CHECK(fil!=NULL); } else if(!strcmp(word,"filament")) { // filament @@ -1494,12 +1493,13 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { CHECKS(itct==2,"filament format: filament_name statement_name statement_text"); line2=strnword(line2,3); CHECKS(line2,"filament format: filament_name statement_name statement_text"); - f=stringfind(sim->filss->fnames,sim->filss->nfil,nm); - CHECKS(f>=0,"filament is unrecognized"); - fil=sim->filss->fillist[f]; - fil=filreadstring(sim,pfp,fil,nm1,line2); + f=filGetFilIndex(sim,nm,&ft); + CHECKS(!(f==-2),"multiple filaments have the same name"); + CHECKS(f>=0,"filament name is unrecognized"); + filtype=sim->filss->filtypes[ft]; + fil=filtype->fillist[f]; + fil=filreadstring(sim,pfp,fil,filtype,nm1,line2); CHECK(fil!=NULL); } -*/ else if(!strcmp(word,"random_filament")) { // random_filament CHECKS(sim->filss,"need to enter a filament type before random_filament"); @@ -1509,7 +1509,8 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { ft=stringfind(filss->ftnames,filss->ntype,nm1); CHECKS(ft>=0,"filament type is unknown"); filtype=filss->filtypes[ft]; - fil=filAddFilament(filtype,NULL,nm); + CHECKS(filtype->klen==-1 || filtype->klen>0,"cannot compute random segments because the filament type has length force constant equals 0"); + fil=filAddFilament(filtype,nm); CHECKS(fil,"unable to add filament to simulation"); line2=strnword(line2,3); @@ -1518,7 +1519,7 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { CHECKS(itct==1,"random_filament format: number [x y z theta phi chi] [thickness]"); CHECKS(i1>0,"number needs to be >0"); line2=strnword(line2,2); - if(fil->nbs==0) { + if(fil->nseg==0) { CHECKS(line2,"missing position and angle information"); itct=sscanf(line2,"%s %s %s %s %s %s",str1,str2,str3,str4,str5,str6); CHECKS(itct==6,"random_filament format: number [x y z theta phi chi] [thickness]"); @@ -1536,10 +1537,7 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { CHECKS(itct==1,"random_segments format: number [x y z theta phi chi] [thickness]"); CHECKS(thick>0,"thickness needs to be >0"); line2=strnword(line2,2); } - if(filtype->isbead) - er=filAddRandomBeads(fil,i1,str1,str2,str3); - else - er=filAddRandomSegments(fil,i1,str1,str2,str3,str4,str5,str6,thick); + er=filAddRandomSegments(fil,i1,str1,str2,str3,str4,str5,str6,thick); CHECKS(er!=2,"random_filament positions need to be 'u' or value"); CHECKS(er!=3,"random_filament angles need to be 'u' or value"); CHECKS(er==0,"BUG: error in filAddRandomSegments"); @@ -2232,7 +2230,7 @@ int loadsim(simptr sim,const char *fileroot,const char *filename,const char *fla er=filloadtype(sim,&pfp,line2); } else if(!strcmp(word,"start_filament")) { // start_filament - er=filloadfil(sim,&pfp,line2,NULL); } + er=filloadfil(sim,&pfp,line2); } else if(!strcmp(word,"start_rules")) { // start_rules CHECKS(0,"Moleculizer support has been discontinued in Smoldyn"); } @@ -2306,7 +2304,7 @@ int simupdate(simptr sim) { if(sim->condition==SCinit && sim->filss) simLog(sim,2," setting up filaments\n"); - er=filsupdate(sim); + er=filupdate(sim); CHECK(er!=1); if(sim->condition==SCinit && sim->graphss)