Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2020-07-06 384af9d04809481663a9fb350212b5e5880955aa
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;
d7639a 111
384af9 112
SP 113     ts_double Se11, Se21, Se22, Se31, Se32, Se33;
114     ts_double Pv11, Pv21, Pv22, Pv31, Pv32, Pv33;
115     ts_double We;
116     ts_double Av, We_Av;
117
118     gsl_matrix *gsl_Sv=gsl_matrix_alloc(3,3);
119     gsl_vector_complex *Sv_eigen=gsl_vector_complex_alloc(3);
120     gsl_eigen_nonsymm_workspace *workspace=gsl_eigen_nonsymm_alloc(3);
121
122     ts_double mprod[7], phi[7], he[7];
123     ts_double Sv[3][3]={{0,0,0},{0,0,0},{0,0,0}};
7d84ef 124     // Here edge vector is calculated
SP 125 //    fprintf(stderr, "Vertex has neighbours=%d\n", vtx->neigh_no);
126     for(jj=0;jj<vtx->neigh_no;jj++){
127     edge_vector_x[jj]=vtx->neigh[jj]->x-vtx->x;
128     edge_vector_y[jj]=vtx->neigh[jj]->y-vtx->y;
129     edge_vector_z[jj]=vtx->neigh[jj]->z-vtx->z;
384af9 130     Av=0;
SP 131     for(i=0; i<vtx->tristar_no; i++){
132         vertex_normal_x=vertex_normal_x + vtx->tristar[i]->xnorm*vtx->tristar[i]->area;
133         vertex_normal_y=vertex_normal_y + vtx->tristar[i]->ynorm*vtx->tristar[i]->area;
134         vertex_normal_z=vertex_normal_z + vtx->tristar[i]->znorm*vtx->tristar[i]->area;
135         Av+=vtx->tristar[i]->area/3;
136     }
137
138     Pv11=1-vertex_normal_x*vertex_normal_x;
139     Pv22=1-vertex_normal_y*vertex_normal_y;
140     Pv33=1-vertex_normal_z*vertex_normal_z;
141     Pv21=vertex_normal_x*vertex_normal_y;
142     Pv31=vertex_normal_x*vertex_normal_z;
143     Pv32=vertex_normal_y*vertex_normal_z;
144
145 //    printf("(%f %f %f)\n", vertex_normal_x, vertex_normal_y, vertex_normal_z);
608cbe 146
SP 147
148     it=vtx;
149     k=vtx->neigh[jj];
150     nei=0;
151         for(i=0;i<it->neigh_no;i++){ // Finds the nn of it, that is k 
152             if(it->neigh[i]==k){
153                 nei=i;
154                 break;
155             }
156         }
157         neip=nei+1;  // I don't like it.. Smells like I must have it in correct order
158         neim=nei-1;
159         if(neip>=it->neigh_no) neip=0;
160         if((ts_int)neim<0) neim=it->neigh_no-1; /* casting is essential... If not
161 there the neim is never <0 !!! */
162   //  fprintf(stderr,"The numbers are: %u %u\n",neip, neim);
163         km=it->neigh[neim];  // We located km and kp
164         kp=it->neigh[neip];
165
166         if(km==NULL || kp==NULL){
384af9 167             fatal("energy_vertex: cannot determine km and kp!",233);
608cbe 168         }
SP 169
170    for(i=0;i<it->tristar_no;i++){
171         for(j=0;j<k->tristar_no;j++){
172             if((it->tristar[i] == k->tristar[j])){ //ce gre za skupen trikotnik
173                 if((it->tristar[i]->vertex[0] == km || it->tristar[i]->vertex[1]
174 == km || it->tristar[i]->vertex[2]== km )){
175                 lm=it->tristar[i];
176          //       lmidx=i;
177                 }
178                 else
179                 {
180                 lp=it->tristar[i];
181          //       lpidx=i;
182                 }
183
184             }
185         }
186     }
384af9 187 if(lm==NULL || lp==NULL) fatal("energy_vertex: Cannot find triangles lm and lp!",233);
d7639a 188
608cbe 189     sumnorm=sqrt( pow((lm->xnorm + lp->xnorm),2) + pow((lm->ynorm + lp->ynorm), 2) + pow((lm->znorm + lp->znorm), 2));
SP 190
191     edge_normal_x[jj]=(lm->xnorm+ lp->xnorm)/sumnorm;
192     edge_normal_y[jj]=(lm->ynorm+ lp->ynorm)/sumnorm;
193     edge_normal_z[jj]=(lm->znorm+ lp->znorm)/sumnorm;
7d84ef 194
SP 195
196     edge_binormal_x[jj]=(edge_normal_y[jj]*edge_vector_z[jj])-(edge_normal_z[jj]*edge_vector_y[jj]);
197     edge_binormal_y[jj]=-(edge_normal_x[jj]*edge_vector_z[jj])+(edge_normal_z[jj]*edge_vector_x[jj]);
198     edge_binormal_z[jj]=(edge_normal_x[jj]*edge_vector_y[jj])-(edge_normal_y[jj]*edge_vector_x[jj]);
199
384af9 200
SP 201     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]);
202     phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm),mprod[jj])+M_PI;
203     he[jj]=2.0*sqrt( pow((edge_vector_x[jj]*2),2) + pow((edge_vector_y[jj]*2), 2) + pow((edge_vector_z[jj]*2), 2))*cos(phi[jj]/2.0);
204
205
206     Se11=edge_binormal_x[jj]*edge_binormal_x[jj]*he[jj];
207     Se21=edge_binormal_x[jj]*edge_binormal_y[jj]*he[jj];
208     Se22=edge_binormal_y[jj]*edge_binormal_y[jj]*he[jj];
209     Se31=edge_binormal_x[jj]*edge_binormal_z[jj]*he[jj];
210     Se32=edge_binormal_y[jj]*edge_binormal_z[jj]*he[jj];
211     Se33=edge_binormal_z[jj]*edge_binormal_z[jj]*he[jj];
212
213     We=vertex_normal_x*edge_normal_x[jj]+vertex_normal_y*edge_normal_y[jj]+vertex_normal_z*edge_normal_z[jj];
214     We_Av=We/Av;
215
216     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) );
217     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));
218     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));
219     
220     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));
221     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));
222     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));
223
224     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));
225     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));
226     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));
227 //    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 228
SP 229     }
384af9 230
SP 231     gsl_matrix_set(gsl_Sv, 0,0, Sv[0][0]);
232     gsl_matrix_set(gsl_Sv, 0,1, Sv[0][1]);
233     gsl_matrix_set(gsl_Sv, 0,2, Sv[0][2]);
234     gsl_matrix_set(gsl_Sv, 1,0, Sv[1][0]);
235     gsl_matrix_set(gsl_Sv, 1,1, Sv[1][1]);
236     gsl_matrix_set(gsl_Sv, 1,2, Sv[1][2]);
237     gsl_matrix_set(gsl_Sv, 2,0, Sv[2][0]);
238     gsl_matrix_set(gsl_Sv, 2,1, Sv[2][1]);
239     gsl_matrix_set(gsl_Sv, 2,2, Sv[2][2]);
240
241     gsl_eigen_nonsymm_params(0, 1, workspace);
242     gsl_eigen_nonsymm(gsl_Sv, Sv_eigen, workspace);
243
244     printf("Eigenvalues: %f+ i%f, %f+i%f, %f+i%f\n", 
245     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 0)), GSL_IMAG(gsl_vector_complex_get(Sv_eigen, 0)),
246     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 1)), GSL_IMAG(gsl_vector_complex_get(Sv_eigen, 1)),
247     GSL_REAL(gsl_vector_complex_get(Sv_eigen, 2)), GSL_IMAG(gsl_vector_complex_get(Sv_eigen, 2))
248     );
7d84ef 249     vtx->energy=0.0;
384af9 250
SP 251     gsl_matrix_free(gsl_Sv);
252     gsl_vector_complex_free(Sv_eigen);
253     gsl_eigen_nonsymm_free(workspace);
7d84ef 254     return TS_SUCCESS;
d7639a 255 }
e5858f 256
SP 257 ts_bool sweep_attraction_bond_energy(ts_vesicle *vesicle){
258     int i;
259     for(i=0;i<vesicle->blist->n;i++){
260         attraction_bond_energy(vesicle->blist->bond[i], vesicle->tape->w);
261     }
262     return TS_SUCCESS;
263 }
264
265
266 inline ts_bool attraction_bond_energy(ts_bond *bond, ts_double w){
267
268     if(fabs(bond->vtx1->c)>1e-16 && fabs(bond->vtx2->c)>1e-16){
032273 269         bond->energy=-w;
e5858f 270     }
SP 271     else {
272         bond->energy=0.0;
273     }
274     return TS_SUCCESS;
275 }
250de4 276
SP 277 ts_double direct_force_energy(ts_vesicle *vesicle, ts_vertex *vtx, ts_vertex *vtx_old){
278     if(fabs(vtx->c)<1e-15) return 0.0;
279 //    printf("was here");
280     if(fabs(vesicle->tape->F)<1e-15) return 0.0;
281
282     ts_double norml,ddp=0.0;
283     ts_uint i;
284     ts_double xnorm=0.0,ynorm=0.0,znorm=0.0;
02d65c 285     /*find normal of the vertex as sum of all the normals of the triangles surrounding it. */
250de4 286     for(i=0;i<vtx->tristar_no;i++){
02d65c 287             xnorm+=vtx->tristar[i]->xnorm;
MF 288             ynorm+=vtx->tristar[i]->ynorm;
289             znorm+=vtx->tristar[i]->znorm;
250de4 290     }
SP 291     /*normalize*/
292     norml=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
293     xnorm/=norml;
294     ynorm/=norml;
295     znorm/=norml;
296     /*calculate ddp, perpendicular displacement*/
c372c1 297     ddp=xnorm*(vtx->x-vtx_old->x)+ynorm*(vtx->y-vtx_old->y)+znorm*(vtx->z-vtx_old->z);
250de4 298     /*calculate dE*/
SP 299 //    printf("ddp=%e",ddp);
300     return vesicle->tape->F*ddp;        
301     
302 }
7ec6fb 303
SP 304 void stretchenergy(ts_vesicle *vesicle, ts_triangle *triangle){
04694f 305     triangle->energy=vesicle->tape->xkA0/2.0*pow((triangle->area/vesicle->tlist->a0-1.0),2);
7ec6fb 306 }