| 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){ |
| 7d84ef | 91 | ts_uint jj, i, j, cnt=0; |
| SP | 92 | ts_double edge_vector_x[7]={0,0,0,0,0,0,0}; |
| 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; | |
| 104 | ts_triangle *triedge[2]={NULL,NULL}; | |
| a63f17 | 105 | |
| 7d84ef | 106 | ts_double sumnorm; |
| d7639a | 107 | |
| 7d84ef | 108 | // Here edge vector is calculated |
| SP | 109 | // fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no); |
| 110 | for(jj=0;jj<vtx->neigh_no;jj++){ | |
| 111 | edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x; | |
| 112 | edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y; | |
| 113 | edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z; | |
| 114 | // We find lm and lp from k->tristar ! | |
| 115 | cnt=0; | |
| 116 | for(i=0;i<vtx->tristar_no;i++){ | |
| 117 | for(j=0;j<vtx->neigh[jj]->tristar_no;j++){ | |
| 118 | if((vtx->tristar[i] == vtx->neigh[jj]->tristar[j])){ //ce gre za skupen trikotnik | |
| 119 | triedge[cnt]=vtx->tristar[i]; | |
| 120 | cnt++; | |
| 121 | } | |
| 122 | } | |
| 123 | } | |
| 124 | if(cnt!=2) fatal("ts_energy_vertex: both triangles not found!", 133); | |
| 125 | sumnorm=sqrt( pow((triedge[0]->xnorm + triedge[1]->xnorm),2) + pow((triedge[0]->ynorm + triedge[1]->ynorm), 2) + pow((triedge[0]->znorm + triedge[1]->znorm), 2)); | |
| d7639a | 126 | |
| 7d84ef | 127 | edge_normal_x[jj]=(triedge[0]->xnorm+ triedge[1]->xnorm)/sumnorm; |
| SP | 128 | edge_normal_y[jj]=(triedge[0]->ynorm+ triedge[1]->ynorm)/sumnorm; |
| 129 | edge_normal_z[jj]=(triedge[0]->znorm+ triedge[1]->znorm)/sumnorm; | |
| 130 | ||
| 131 | ||
| 132 | edge_binormal_x[jj]=(edge_normal_y[jj]*edge_vector_z[jj])-(edge_normal_z[jj]*edge_vector_y[jj]); | |
| 133 | edge_binormal_y[jj]=-(edge_normal_x[jj]*edge_vector_z[jj])+(edge_normal_z[jj]*edge_vector_x[jj]); | |
| 134 | edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]); | |
| 135 | ||
| 136 | printf("(%f %f %f); (%f %f %f); (%f %f %f), %d\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],triedge[0]->idx); | |
| 137 | ||
| 138 | } | |
| 139 | for(i=0; i<vtx->tristar_no; i++){ | |
| 140 | vertex_normal_x=vertex_normal_x + vtx->tristar[i]->xnorm*vtx->tristar[i]->area; | |
| 141 | vertex_normal_y=vertex_normal_y + vtx->tristar[i]->ynorm*vtx->tristar[i]->area; | |
| 142 | vertex_normal_z=vertex_normal_z + vtx->tristar[i]->znorm*vtx->tristar[i]->area; | |
| 143 | } | |
| 144 | printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z); | |
| 145 | vtx->energy=0.0; | |
| 146 | return TS_SUCCESS; | |
| d7639a | 147 | } |
| e5858f | 148 | |
| SP | 149 | ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){ |
| 150 | int i; | |
| 151 | for(i=0;i<vesicle->blist->n;i++){ | |
| 152 | attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w); | |
| 153 | } | |
| 154 | return TS_SUCCESS; | |
| 155 | } | |
| 156 | ||
| 157 | ||
| 158 | inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){ | |
| 159 | ||
| 160 | if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){ | |
| 032273 | 161 | bond->energy=-w; |
| e5858f | 162 | } |
| SP | 163 | else { |
| 164 | bond->energy=0.0; | |
| 165 | } | |
| 166 | return TS_SUCCESS; | |
| 167 | } | |
| 250de4 | 168 | |
| SP | 169 | ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){ |
| 170 | if(fabs(vtx->c)<1e-15) return 0.0; | |
| 171 | // printf("was here"); | |
| 172 | if(fabs(vesicle->tape->F)<1e-15) return 0.0; | |
| 173 | ||
| 174 | ts_double norml,ddp=0.0; | |
| 175 | ts_uint i; | |
| 176 | ts_double xnorm=0.0,ynorm=0.0,znorm=0.0; | |
| 02d65c | 177 | /*find normal of the vertex as sum of all the normals of the triangles surrounding it. */ |
| 250de4 | 178 | for(i=0;i<vtx->tristar_no;i++){ |
| 02d65c | 179 | xnorm+=vtx->tristar[i]->xnorm; |
| MF | 180 | ynorm+=vtx->tristar[i]->ynorm; |
| 181 | znorm+=vtx->tristar[i]->znorm; | |
| 250de4 | 182 | } |
| SP | 183 | /*normalize*/ |
| 184 | norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm); | |
| 185 | xnorm/=norml; | |
| 186 | ynorm/=norml; | |
| 187 | znorm/=norml; | |
| 188 | /*calculate ddp, perpendicular displacement*/ | |
| c372c1 | 189 | ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z); |
| 250de4 | 190 | /*calculate dE*/ |
| SP | 191 | // printf("ddp=%e",ddp); |
| 192 | return vesicle->tape->F*ddp; | |
| 193 | ||
| 194 | } | |
| 7ec6fb | 195 | |
| SP | 196 | void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){ |
| 04694f | 197 | triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2); |
| 7ec6fb | 198 | } |
