Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2020-10-07 8493be19b22681868f4f6686da592f4d149cdf05
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>
384af9 9 #include <gsl/gsl_vector_complex.h>
SP 10 #include <gsl/gsl_matrix.h>
11 #include <gsl/gsl_eigen.h>
a00f10 12 /** @brief Wrapper that calculates energy of every vertex in vesicle
SP 13  *  
14  *  Function calculated energy of every vertex in vesicle. It can be used in
15  *  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. 
16  *  @param *vesicle is a pointer to vesicle.
17  *  @returns TS_SUCCESS on success.
18 */
d7639a 19 ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){
SP 20
f74313 21     ts_uint i;
d7639a 22     
f74313 23     ts_vertex_list *vlist=vesicle->vlist;
SP 24     ts_vertex **vtx=vlist->vtx;
d7639a 25
SP 26     for(i=0;i<vlist->n;i++){
f74313 27         energy_vertex(vtx[i]);
b01cc1 28         
d7639a 29     }
SP 30
31     return TS_SUCCESS;
32 }
33
a00f10 34 /** @brief Calculate energy of a bond (in models where energy is bond related)
SP 35  *
36  *  This function is experimental and currently only used in polymeres calculation (PEGs or polymeres inside the vesicle).
37  *
38  *  @param *bond is a pointer to a bond between two vertices in polymere
39  *  @param *poly is a pointer to polymere in which we calculate te energy of the bond
40  *  @returns TS_SUCCESS on successful calculation
41 */
fedf2b 42 inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){
304510 43 //TODO: This value to be changed and implemented in data structure:
M 44     ts_double d_relaxed=1.0;
45     bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2);
fedf2b 46     return TS_SUCCESS;
M 47 };
48
e6efc6 49 /** @brief Calculation of the bending energy of the vertex.
a00f10 50  *  
e6efc6 51  *  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 52  * \f$c/2\f$ is spontaneous curvature and \f$s\f$ is area per vertex \f$i\f$.
53  *
54  * Nearest neighbors (NN) must be ordered in counterclockwise direction for this function to work.
a00f10 55  *  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 56  *  
854cb6 57 \begin{tikzpicture}{
a00f10 58 \coordinate[label=below:$i$] (i) at (2,0);
SP 59 \coordinate[label=left:$j_m$] (jm) at (0,3.7);
60 \coordinate[label=above:$j$] (j) at (2.5,6.4);
61 \coordinate[label=right:$j_p$] (jp) at (4,2.7);
d7639a 62
a00f10 63 \draw (i) -- (jm) -- (j) -- (jp) -- (i) -- (j);
SP 64
65 \begin{scope}
66 \path[clip] (jm)--(i)--(j);
67 \draw (jm) circle (0.8);
68 \node[right] at (jm) {$\varphi_m$};
69 \end{scope}
70
71 \begin{scope}
72 \path[clip] (jp)--(i)--(j);
73 \draw (jp) circle (0.8);
74 \node[left] at (jp) {$\varphi_p$};
75 \end{scope}
76
77 %%vertices
78 \draw [fill=gray] (i) circle (0.1);
79 \draw [fill=white] (j) circle (0.1);
80 \draw [fill=white] (jp) circle (0.1);
81 \draw [fill=white] (jm) circle (0.1);
82 %\node[draw,circle,fill=white] at (i) {};
854cb6 83 \end{tikzpicture}
a00f10 84
SP 85  * 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).
86  *
87  * From the curvature the enery is calculated by equation \f$E=\frac{1}{2}\mathbf{h}\cdot\mathbf{h}\f$.
88  * @param *vtx is a pointer to vertex at which we want to calculate the energy
89  * @returns TS_SUCCESS on successful calculation.
90 */
d7639a 91 inline ts_bool energy_vertex(ts_vertex *vtx){
608cbe 92     ts_uint jj, i, j;
7d84ef 93     ts_double edge_vector_x[7]={0,0,0,0,0,0,0};
SP 94     ts_double edge_vector_y[7]={0,0,0,0,0,0,0};
95     ts_double edge_vector_z[7]={0,0,0,0,0,0,0};
96     ts_double edge_normal_x[7]={0,0,0,0,0,0,0};
97     ts_double edge_normal_y[7]={0,0,0,0,0,0,0};
98     ts_double edge_normal_z[7]={0,0,0,0,0,0,0};
99     ts_double edge_binormal_x[7]={0,0,0,0,0,0,0};
100     ts_double edge_binormal_y[7]={0,0,0,0,0,0,0};
101     ts_double edge_binormal_z[7]={0,0,0,0,0,0,0};
102     ts_double vertex_normal_x=0.0;
103     ts_double vertex_normal_y=0.0;
104     ts_double vertex_normal_z=0.0;
608cbe 105 //    ts_triangle *triedge[2]={NULL,NULL};
a63f17 106
608cbe 107     ts_uint nei,neip,neim;
SP 108     ts_vertex *it, *k, *kp,*km;
109     ts_triangle *lm=NULL, *lp=NULL;
7d84ef 110     ts_double sumnorm;
8493be 111     ts_double temp_length;
d7639a 112
384af9 113
SP 114     ts_double Se11, Se21, Se22, Se31, Se32, Se33;
115     ts_double Pv11, Pv21, Pv22, Pv31, Pv32, Pv33;
116     ts_double We;
117     ts_double Av, We_Av;
118
119     gsl_matrix *gsl_Sv=gsl_matrix_alloc(3,3);
120     gsl_vector_complex *Sv_eigen=gsl_vector_complex_alloc(3);
121     gsl_eigen_nonsymm_workspace *workspace=gsl_eigen_nonsymm_alloc(3);
122
123     ts_double mprod[7], phi[7], he[7];
124     ts_double Sv[3][3]={{0,0,0},{0,0,0},{0,0,0}};
7d84ef 125     // Here edge vector is calculated
SP 126 //    fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no);
8493be 127
SP 128
129
130
384af9 131     Av=0;
SP 132     for(i=0; i<vtx->tristar_no; i++){
133         vertex_normal_x=vertex_normal_x + vtx->tristar[i]->xnorm*vtx->tristar[i]->area;
134         vertex_normal_y=vertex_normal_y + vtx->tristar[i]->ynorm*vtx->tristar[i]->area;
135         vertex_normal_z=vertex_normal_z + vtx->tristar[i]->znorm*vtx->tristar[i]->area;
136         Av+=vtx->tristar[i]->area/3;
137     }
8493be 138     temp_length=sqrt(pow(vertex_normal_x,2)+pow(vertex_normal_y,2)+pow(vertex_normal_z,2));
SP 139     vertex_normal_x=vertex_normal_x/temp_length;
140     vertex_normal_y=vertex_normal_y/temp_length;
141     vertex_normal_z=vertex_normal_z/temp_length;
384af9 142
SP 143     Pv11=1-vertex_normal_x*vertex_normal_x;
144     Pv22=1-vertex_normal_y*vertex_normal_y;
145     Pv33=1-vertex_normal_z*vertex_normal_z;
146     Pv21=vertex_normal_x*vertex_normal_y;
147     Pv31=vertex_normal_x*vertex_normal_z;
148     Pv32=vertex_normal_y*vertex_normal_z;
149
8493be 150
SP 151
152
153     for(jj=0;jj<vtx->neigh_no;jj++){
154     edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x;
155     edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y;
156     edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z;
157
158     //Here we calculate normalized edge vector
159
160     temp_length=sqrt(edge_vector_x[jj]*edge_vector_x[jj]+edge_vector_y[jj]*edge_vector_y[jj]+edge_vector_z[jj]*edge_vector_z[jj]);
161     edge_vector_x[jj]=edge_vector_x[jj]/temp_length;
162     edge_vector_y[jj]=edge_vector_y[jj]/temp_length;
163     edge_vector_z[jj]=edge_vector_z[jj]/temp_length;
164
165     //end normalization
384af9 166 //    printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z);
608cbe 167
SP 168
169     it=vtx;
170     k=vtx->neigh[jj];
171     nei=0;
172         for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k 
173             if(it->neigh[i]==k){
174                 nei=i;
175                 break;
176             }
177         }
178         neip=nei+1;  // I don't like it.. Smells like I must have it in correct order
179         neim=nei-1;
180         if(neip>=it->neigh_no) neip=0;
181         if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not
182 there the neim is never <0 !!! */
183   //  fprintf(stderr,"The numbers are: %u %u\n",neip, neim);
184         km=it->neigh[neim];  // We located km and kp
185         kp=it->neigh[neip];
186
187         if(km==NULL || kp==NULL){
384af9 188             fatal("energy_vertex: cannot determine km and kp!",233);
608cbe 189         }
SP 190
191    for(i=0;i<it->tristar_no;i++){
192         for(j=0;j<k->tristar_no;j++){
193             if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
194                 if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
195 == km || it->tristar[i]->vertex[2]== km )){
196                 lm=it->tristar[i];
197          //       lmidx=i;
198                 }
199                 else
200                 {
201                 lp=it->tristar[i];
202          //       lpidx=i;
203                 }
204
205             }
206         }
207     }
384af9 208 if(lm==NULL || lp==NULL) fatal("energy_vertex: Cannot find triangles lm and lp!",233);
d7639a 209
8493be 210     //Triangle normals are NORMALIZED!
SP 211
608cbe 212     sumnorm=sqrt( pow((lm->xnorm + lp->xnorm),2) + pow((lm->ynorm + lp->ynorm), 2) + pow((lm->znorm + lp->znorm), 2));
SP 213
214     edge_normal_x[jj]=(lm->xnorm+ lp->xnorm)/sumnorm;
215     edge_normal_y[jj]=(lm->ynorm+ lp->ynorm)/sumnorm;
216     edge_normal_z[jj]=(lm->znorm+ lp->znorm)/sumnorm;
7d84ef 217
SP 218
219     edge_binormal_x[jj]=(edge_normal_y[jj]*edge_vector_z[jj])-(edge_normal_z[jj]*edge_vector_y[jj]);
220     edge_binormal_y[jj]=-(edge_normal_x[jj]*edge_vector_z[jj])+(edge_normal_z[jj]*edge_vector_x[jj]);
221     edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]);
222
384af9 223
SP 224     mprod[jj]=it->x*(k->y*edge_vector_z[jj]-edge_vector_y[jj]*k->z)-it->y*(k->x*edge_vector_z[jj]-k->z*edge_vector_x[jj])+it->z*(k->x*edge_vector_y[jj]-k->y*edge_vector_x[jj]);
8493be 225     phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-8),mprod[jj])+M_PI;
SP 226 //    printf("ACOS arg=%e\n", lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm);
227     //he was multiplied with 2 before...
228     he[jj]=sqrt( pow((edge_vector_x[jj]),2) + pow((edge_vector_y[jj]), 2) + pow((edge_vector_z[jj]), 2))*cos(phi[jj]/2.0);
229 //    printf("phi[%d]=%f\n", jj,phi[jj]);
384af9 230
SP 231     Se11=edge_binormal_x[jj]*edge_binormal_x[jj]*he[jj];
232     Se21=edge_binormal_x[jj]*edge_binormal_y[jj]*he[jj];
233     Se22=edge_binormal_y[jj]*edge_binormal_y[jj]*he[jj];
234     Se31=edge_binormal_x[jj]*edge_binormal_z[jj]*he[jj];
235     Se32=edge_binormal_y[jj]*edge_binormal_z[jj]*he[jj];
236     Se33=edge_binormal_z[jj]*edge_binormal_z[jj]*he[jj];
237
238     We=vertex_normal_x*edge_normal_x[jj]+vertex_normal_y*edge_normal_y[jj]+vertex_normal_z*edge_normal_z[jj];
239     We_Av=We/Av;
240
241     Sv[0][0]+=We_Av* ( Pv11*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv21*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv31*(Pv11*Se31+Pv21*Se32+Pv31*Se33) );
242     Sv[0][1]+=We_Av* (Pv21*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv22*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv32*(Pv11*Se31+Pv21*Se32+Pv31*Se33));
243     Sv[0][2]+=We_Av* (Pv31*(Pv11*Se11+Pv21*Se21+Pv31*Se31)+Pv32*(Pv11*Se21+Pv21*Se22+Pv31*Se32)+Pv33*(Pv11*Se31+Pv21*Se32+Pv31*Se33));
244     
245     Sv[1][0]+=We_Av* (Pv11*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv21*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv31*(Pv21*Se31+Pv22*Se32+Pv32*Se33));
246     Sv[1][1]+=We_Av* (Pv21*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv22*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv32*(Pv21*Se31+Pv22*Se32+Pv32*Se33));
247     Sv[1][2]+=We_Av* (Pv31*(Pv21*Se11+Pv22*Se21+Pv32*Se31)+Pv32*(Pv21*Se21+Pv22*Se22+Pv32*Se32)+Pv33*(Pv21*Se31+Pv22*Se32+Pv32*Se33));
248
249     Sv[2][0]+=We_Av* (Pv11*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv21*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv31*(Pv31*Se31+Pv32*Se32+Pv33*Se33));
250     Sv[2][1]+=We_Av* (Pv21*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv22*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv32*(Pv31*Se31+Pv32*Se32+Pv33*Se33));
251     Sv[2][2]+=We_Av* (Pv31*(Pv31*Se11+Pv32*Se21+Pv33*Se31)+Pv32*(Pv31*Se21+Pv32*Se22+Pv33*Se32)+Pv33*(Pv31*Se31+Pv32*Se32+Pv33*Se33));
252 //    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 253
8493be 254     } // END FOR JJ
384af9 255
SP 256     gsl_matrix_set(gsl_Sv, 0,0, Sv[0][0]);
257     gsl_matrix_set(gsl_Sv, 0,1, Sv[0][1]);
258     gsl_matrix_set(gsl_Sv, 0,2, Sv[0][2]);
259     gsl_matrix_set(gsl_Sv, 1,0, Sv[1][0]);
260     gsl_matrix_set(gsl_Sv, 1,1, Sv[1][1]);
261     gsl_matrix_set(gsl_Sv, 1,2, Sv[1][2]);
262     gsl_matrix_set(gsl_Sv, 2,0, Sv[2][0]);
263     gsl_matrix_set(gsl_Sv, 2,1, Sv[2][1]);
264     gsl_matrix_set(gsl_Sv, 2,2, Sv[2][2]);
265
8493be 266 //    printf("Se= %f, %f, %f\n    %f, %f, %f\n    %f, %f, %f\n", Se11, Se21, Se31, Se21, Se22, Se32, Se31, Se32, Se33);
SP 267 //    printf("Pv= %f, %f, %f\n    %f, %f, %f\n    %f, %f, %f\n", Pv11, Pv21, Pv31, Pv21, Pv22, Pv32, Pv31, Pv32, Pv33);
268     printf("Sv= %f, %f, %f\n    %f, %f, %f\n    %f, %f, %f\n", Sv[0][0], Sv[0][1], Sv[0][2], Sv[1][0], Sv[1][1], Sv[1][2], Sv[2][0], Sv[2][1], Sv[2][2]);
269
384af9 270     gsl_eigen_nonsymm_params(0, 1, workspace);
SP 271     gsl_eigen_nonsymm(gsl_Sv, Sv_eigen, workspace);
272
8493be 273     printf("Eigenvalues: %f, %f, %f\n", 
SP 274     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 0)),
275     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 1)),
276     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 2))
384af9 277     );
7d84ef 278     vtx->energy=0.0;
384af9 279
SP 280     gsl_matrix_free(gsl_Sv);
281     gsl_vector_complex_free(Sv_eigen);
282     gsl_eigen_nonsymm_free(workspace);
7d84ef 283     return TS_SUCCESS;
d7639a 284 }
e5858f 285
SP 286 ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){
287     int i;
288     for(i=0;i<vesicle->blist->n;i++){
289         attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w);
290     }
291     return TS_SUCCESS;
292 }
293
294
295 inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){
296
297     if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){
032273 298         bond->energy=-w;
e5858f 299     }
SP 300     else {
301         bond->energy=0.0;
302     }
303     return TS_SUCCESS;
304 }
250de4 305
SP 306 ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){
307     if(fabs(vtx->c)<1e-15) return 0.0;
308 //    printf("was here");
309     if(fabs(vesicle->tape->F)<1e-15) return 0.0;
310
311     ts_double norml,ddp=0.0;
312     ts_uint i;
313     ts_double xnorm=0.0,ynorm=0.0,znorm=0.0;
02d65c 314     /*find normal of the vertex as sum of all the normals of the triangles surrounding it. */
250de4 315     for(i=0;i<vtx->tristar_no;i++){
02d65c 316             xnorm+=vtx->tristar[i]->xnorm;
MF 317             ynorm+=vtx->tristar[i]->ynorm;
318             znorm+=vtx->tristar[i]->znorm;
250de4 319     }
SP 320     /*normalize*/
321     norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
322     xnorm/=norml;
323     ynorm/=norml;
324     znorm/=norml;
325     /*calculate ddp, perpendicular displacement*/
c372c1 326     ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z);
250de4 327     /*calculate dE*/
SP 328 //    printf("ddp=%e",ddp);
329     return vesicle->tape->F*ddp;        
330     
331 }
7ec6fb 332
SP 333 void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){
04694f 334     triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2);
7ec6fb 335 }