| commit | author | age | ||
| 7f6076 | 1 | /* vim: set ts=4 sts=4 sw=4 noet : */ |
| d7639a | 2 | #include<stdlib.h> |
| SP | 3 | #include "general.h" |
| 4 | #include "energy.h" | |
| 5 | #include "vertex.h" | |
| 7d84ef | 6 | #include "bond.h" |
| d7639a | 7 | #include<math.h> |
| SP | 8 | #include<stdio.h> |
| a00f10 | 9 | |
| SP | 10 | |
| 11 | /** @brief Wrapper that calculates energy of every vertex in vesicle | |
| 12 | * | |
| 13 | * Function calculated energy of every vertex in vesicle. It can be used in | |
| 14 | * initialization procedure or in recalculation of the energy after non-MCsweep * operations. However, when random move of vertex or flip of random bond occur * call to this function is not necessary nor recommended. | |
| 15 | * @param *vesicle is a pointer to vesicle. | |
| 16 | * @returns TS_SUCCESS on success. | |
| 17 | */ | |
| d7639a | 18 | ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){ |
| SP | 19 | |
| f74313 | 20 | ts_uint i; |
| d7639a | 21 | |
| f74313 | 22 | ts_vertex_list *vlist=vesicle->vlist; |
| SP | 23 | ts_vertex **vtx=vlist->vtx; |
| d7639a | 24 | |
| SP | 25 | for(i=0;i<vlist->n;i++){ |
| f74313 | 26 | energy_vertex(vtx[i]); |
| b01cc1 | 27 | |
| d7639a | 28 | } |
| SP | 29 | |
| 30 | return TS_SUCCESS; | |
| 31 | } | |
| 32 | ||
| a00f10 | 33 | /** @brief Calculate energy of a bond (in models where energy is bond related) |
| SP | 34 | * |
| 35 | * This function is experimental and currently only used in polymeres calculation (PEGs or polymeres inside the vesicle). | |
| 36 | * | |
| 37 | * @param *bond is a pointer to a bond between two vertices in polymere | |
| 38 | * @param *poly is a pointer to polymere in which we calculate te energy of the bond | |
| 39 | * @returns TS_SUCCESS on successful calculation | |
| 40 | */ | |
| fedf2b | 41 | inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){ |
| 304510 | 42 | //TODO: This value to be changed and implemented in data structure: |
| M | 43 | ts_double d_relaxed=1.0; |
| 44 | bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2); | |
| fedf2b | 45 | return TS_SUCCESS; |
| M | 46 | }; |
| 47 | ||
| e6efc6 | 48 | /** @brief Calculation of the bending energy of the vertex. |
| a00f10 | 49 | * |
| e6efc6 | 50 | * Main function that calculates energy of the vertex \f$i\f$. Function returns \f$\frac{1}{2}(c_1+c_2-c)^2 s\f$, where \f$(c_1+c_2)/2\f$ is mean curvature, |
| SP | 51 | * \f$c/2\f$ is spontaneous curvature and \f$s\f$ is area per vertex \f$i\f$. |
| 52 | * | |
| 53 | * Nearest neighbors (NN) must be ordered in counterclockwise direction for this function to work. | |
| a00f10 | 54 | * Firstly NNs that form two neighboring triangles are found (\f$j_m\f$, \f$j_p\f$ and common \f$j\f$). Later, the scalar product of vectors \f$x_1=(\mathbf{i}-\mathbf{j_p})\cdot (\mathbf{i}-\mathbf{j_p})(\mathbf{i}-\mathbf{j_p})\f$, \f$x_2=(\mathbf{j}-\mathbf{j_p})\cdot (\mathbf{j}-\mathbf{j_p})\f$ and \f$x_3=(\mathbf{j}-\mathbf{j_p})\cdot (\mathbf{i}-\mathbf{j_p})\f$ are calculated. From these three vectors the \f$c_{tp}=\frac{1}{\tan(\varphi_p)}\f$ is calculated, where \f$\varphi_p\f$ is the inner angle at vertex \f$j_p\f$. The procedure is repeated for \f$j_m\f$ instead of \f$j_p\f$ resulting in \f$c_{tn}\f$. |
| SP | 55 | * |
| 854cb6 | 56 | \begin{tikzpicture}{ |
| a00f10 | 57 | \coordinate[label=below:$i$] (i) at (2,0); |
| SP | 58 | \coordinate[label=left:$j_m$] (jm) at (0,3.7); |
| 59 | \coordinate[label=above:$j$] (j) at (2.5,6.4); | |
| 60 | \coordinate[label=right:$j_p$] (jp) at (4,2.7); | |
| d7639a | 61 | |
| a00f10 | 62 | \draw (i) -- (jm) -- (j) -- (jp) -- (i) -- (j); |
| SP | 63 | |
| 64 | \begin{scope} | |
| 65 | \path[clip] (jm)--(i)--(j); | |
| 66 | \draw (jm) circle (0.8); | |
| 67 | \node[right] at (jm) {$\varphi_m$}; | |
| 68 | \end{scope} | |
| 69 | ||
| 70 | \begin{scope} | |
| 71 | \path[clip] (jp)--(i)--(j); | |
| 72 | \draw (jp) circle (0.8); | |
| 73 | \node[left] at (jp) {$\varphi_p$}; | |
| 74 | \end{scope} | |
| 75 | ||
| 76 | %%vertices | |
| 77 | \draw [fill=gray] (i) circle (0.1); | |
| 78 | \draw [fill=white] (j) circle (0.1); | |
| 79 | \draw [fill=white] (jp) circle (0.1); | |
| 80 | \draw [fill=white] (jm) circle (0.1); | |
| 81 | %\node[draw,circle,fill=white] at (i) {}; | |
| 854cb6 | 82 | \end{tikzpicture} |
| a00f10 | 83 | |
| SP | 84 | * The curvature is then calculated as \f$\mathbf{h}=\frac{1}{2}\Sigma_{k=0}^{\mathrm{neigh\_no}} c_{tp}^{(k)}+c_{tm}^{(k)} (\mathbf{j_k}-\mathbf{i})\f$, where \f$c_{tp}^{(k)}+c_{tm}^k=2\sigma^{(k)}\f$ (length in dual lattice?) and the previous equation can be written as \f$\mathbf{h}=\Sigma_{k=0}^{\mathrm{neigh\_no}}\sigma^{(k)}\cdot(\mathbf{j}-\mathbf{i})\f$ (See Kroll, p. 384, eq 70). |
| 85 | * | |
| 86 | * From the curvature the enery is calculated by equation \f$E=\frac{1}{2}\mathbf{h}\cdot\mathbf{h}\f$. | |
| 87 | * @param *vtx is a pointer to vertex at which we want to calculate the energy | |
| 88 | * @returns TS_SUCCESS on successful calculation. | |
| 89 | */ | |
| d7639a | 90 | inline ts_bool energy_vertex(ts_vertex *vtx){ |
| 608cbe | 91 | ts_uint jj, i, j; |
| 7d84ef | 92 | ts_double edge_vector_x[7]={0,0,0,0,0,0,0}; |
| SP | 93 | ts_double edge_vector_y[7]={0,0,0,0,0,0,0}; |
| 94 | ts_double edge_vector_z[7]={0,0,0,0,0,0,0}; | |
| 95 | ts_double edge_normal_x[7]={0,0,0,0,0,0,0}; | |
| 96 | ts_double edge_normal_y[7]={0,0,0,0,0,0,0}; | |
| 97 | ts_double edge_normal_z[7]={0,0,0,0,0,0,0}; | |
| 98 | ts_double edge_binormal_x[7]={0,0,0,0,0,0,0}; | |
| 99 | ts_double edge_binormal_y[7]={0,0,0,0,0,0,0}; | |
| 100 | ts_double edge_binormal_z[7]={0,0,0,0,0,0,0}; | |
| 101 | ts_double vertex_normal_x=0.0; | |
| 102 | ts_double vertex_normal_y=0.0; | |
| 103 | ts_double vertex_normal_z=0.0; | |
| 608cbe | 104 | // ts_triangle *triedge[2]={NULL,NULL}; |
| a63f17 | 105 | |
| 608cbe | 106 | ts_uint nei,neip,neim; |
| SP | 107 | ts_vertex *it, *k, *kp,*km; |
| 108 | ts_triangle *lm=NULL, *lp=NULL; | |
| 7d84ef | 109 | ts_double sumnorm; |
| d7639a | 110 | |
| 7d84ef | 111 | // Here edge vector is calculated |
| SP | 112 | // fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no); |
| 113 | for(jj=0;jj<vtx->neigh_no;jj++){ | |
| 114 | edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x; | |
| 115 | edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y; | |
| 116 | edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z; | |
| 608cbe | 117 | |
| SP | 118 | |
| 119 | it=vtx; | |
| 120 | k=vtx->neigh[jj]; | |
| 121 | nei=0; | |
| 122 | for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k | |
| 123 | if(it->neigh[i]==k){ | |
| 124 | nei=i; | |
| 125 | break; | |
| 126 | } | |
| 127 | } | |
| 128 | neip=nei+1; // I don't like it.. Smells like I must have it in correct order | |
| 129 | neim=nei-1; | |
| 130 | if(neip>=it->neigh_no) neip=0; | |
| 131 | if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not | |
| 132 | there the neim is never <0 !!! */ | |
| 133 | // fprintf(stderr,"The numbers are: %u %u\n",neip, neim); | |
| 134 | km=it->neigh[neim]; // We located km and kp | |
| 135 | kp=it->neigh[neip]; | |
| 136 | ||
| 137 | if(km==NULL || kp==NULL){ | |
| 138 | fatal("In bondflip, cannot determine km and kp!",999); | |
| 139 | } | |
| 140 | ||
| 141 | for(i=0;i<it->tristar_no;i++){ | |
| 142 | for(j=0;j<k->tristar_no;j++){ | |
| 143 | if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik | |
| 144 | if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1] | |
| 145 | == km || it->tristar[i]->vertex[2]== km )){ | |
| 146 | lm=it->tristar[i]; | |
| 147 | // lmidx=i; | |
| 148 | } | |
| 149 | else | |
| 150 | { | |
| 151 | lp=it->tristar[i]; | |
| 152 | // lpidx=i; | |
| 153 | } | |
| 154 | ||
| 155 | } | |
| 156 | } | |
| 157 | } | |
| 158 | if(lm==NULL || lp==NULL) fatal("ts_flip_bond: Cannot find triangles lm and lp!",999); | |
| 159 | ||
| 160 | ||
| 161 | /* | |
| 7d84ef | 162 | // We find lm and lp from k->tristar ! |
| SP | 163 | cnt=0; |
| 164 | for(i=0;i<vtx->tristar_no;i++){ | |
| 165 | for(j=0;j<vtx->neigh[jj]->tristar_no;j++){ | |
| 166 | if((vtx->tristar[i] == vtx->neigh[jj]->tristar[j])){ //ce gre za skupen trikotnik | |
| 167 | triedge[cnt]=vtx->tristar[i]; | |
| 168 | cnt++; | |
| 169 | } | |
| 170 | } | |
| 171 | } | |
| 172 | if(cnt!=2) fatal("ts_energy_vertex: both triangles not found!", 133); | |
| 608cbe | 173 | */ |
| d7639a | 174 | |
| 608cbe | 175 | sumnorm=sqrt( pow((lm->xnorm + lp->xnorm),2) + pow((lm->ynorm + lp->ynorm), 2) + pow((lm->znorm + lp->znorm), 2)); |
| SP | 176 | |
| 177 | edge_normal_x[jj]=(lm->xnorm+ lp->xnorm)/sumnorm; | |
| 178 | edge_normal_y[jj]=(lm->ynorm+ lp->ynorm)/sumnorm; | |
| 179 | edge_normal_z[jj]=(lm->znorm+ lp->znorm)/sumnorm; | |
| 7d84ef | 180 | |
| SP | 181 | |
| 182 | edge_binormal_x[jj]=(edge_normal_y[jj]*edge_vector_z[jj])-(edge_normal_z[jj]*edge_vector_y[jj]); | |
| 183 | edge_binormal_y[jj]=-(edge_normal_x[jj]*edge_vector_z[jj])+(edge_normal_z[jj]*edge_vector_x[jj]); | |
| 184 | edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]); | |
| 185 | ||
| 608cbe | 186 | printf("(%f %f %f); (%f %f %f); (%f %f %f)\n", edge_vector_x[jj], edge_vector_y[jj], edge_vector_z[jj], edge_normal_x[jj], edge_normal_y[jj], edge_normal_z[jj], edge_binormal_x[jj], edge_binormal_y[jj], edge_binormal_z[jj]); |
| 7d84ef | 187 | |
| SP | 188 | } |
| 189 | for(i=0; i<vtx->tristar_no; i++){ | |
| 190 | vertex_normal_x=vertex_normal_x + vtx->tristar[i]->xnorm*vtx->tristar[i]->area; | |
| 191 | vertex_normal_y=vertex_normal_y + vtx->tristar[i]->ynorm*vtx->tristar[i]->area; | |
| 192 | vertex_normal_z=vertex_normal_z + vtx->tristar[i]->znorm*vtx->tristar[i]->area; | |
| 193 | } | |
| 194 | printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z); | |
| 195 | vtx->energy=0.0; | |
| 196 | return TS_SUCCESS; | |
| d7639a | 197 | } |
| e5858f | 198 | |
| SP | 199 | ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){ |
| 200 | int i; | |
| 201 | for(i=0;i<vesicle->blist->n;i++){ | |
| 202 | attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w); | |
| 203 | } | |
| 204 | return TS_SUCCESS; | |
| 205 | } | |
| 206 | ||
| 207 | ||
| 208 | inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){ | |
| 209 | ||
| 210 | if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){ | |
| 032273 | 211 | bond->energy=-w; |
| e5858f | 212 | } |
| SP | 213 | else { |
| 214 | bond->energy=0.0; | |
| 215 | } | |
| 216 | return TS_SUCCESS; | |
| 217 | } | |
| 250de4 | 218 | |
| SP | 219 | ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){ |
| 220 | if(fabs(vtx->c)<1e-15) return 0.0; | |
| 221 | // printf("was here"); | |
| 222 | if(fabs(vesicle->tape->F)<1e-15) return 0.0; | |
| 223 | ||
| 224 | ts_double norml,ddp=0.0; | |
| 225 | ts_uint i; | |
| 226 | ts_double xnorm=0.0,ynorm=0.0,znorm=0.0; | |
| 02d65c | 227 | /*find normal of the vertex as sum of all the normals of the triangles surrounding it. */ |
| 250de4 | 228 | for(i=0;i<vtx->tristar_no;i++){ |
| 02d65c | 229 | xnorm+=vtx->tristar[i]->xnorm; |
| MF | 230 | ynorm+=vtx->tristar[i]->ynorm; |
| 231 | znorm+=vtx->tristar[i]->znorm; | |
| 250de4 | 232 | } |
| SP | 233 | /*normalize*/ |
| 234 | norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm); | |
| 235 | xnorm/=norml; | |
| 236 | ynorm/=norml; | |
| 237 | znorm/=norml; | |
| 238 | /*calculate ddp, perpendicular displacement*/ | |
| c372c1 | 239 | ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z); |
| 250de4 | 240 | /*calculate dE*/ |
| SP | 241 | // printf("ddp=%e",ddp); |
| 242 | return vesicle->tape->F*ddp; | |
| 243 | ||
| 244 | } | |
| 7ec6fb | 245 | |
| SP | 246 | void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){ |
| 04694f | 247 | triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2); |
| 7ec6fb | 248 | } |
