Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2016-02-29 bb40336eb8a2acc9f2e3d2e9059b7d90afe5c0c4
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include "general.h"
3 #include "energy.h"
4 #include "vertex.h"
5 #include<math.h>
6 #include<stdio.h>
a00f10 7
SP 8
9 /** @brief Wrapper that calculates energy of every vertex in vesicle
10  *  
11  *  Function calculated energy of every vertex in vesicle. It can be used in
12  *  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. 
13  *  @param *vesicle is a pointer to vesicle.
14  *  @returns TS_SUCCESS on success.
15 */
d7639a 16 ts_bool mean_curvature_and_energy(ts_vesicle *vesicle){
SP 17
f74313 18     ts_uint i;
d7639a 19     
f74313 20     ts_vertex_list *vlist=vesicle->vlist;
SP 21     ts_vertex **vtx=vlist->vtx;
d7639a 22
SP 23     for(i=0;i<vlist->n;i++){
f74313 24         energy_vertex(vtx[i]);
b01cc1 25         
d7639a 26     }
SP 27
28     return TS_SUCCESS;
29 }
30
a00f10 31 /** @brief Calculate energy of a bond (in models where energy is bond related)
SP 32  *
33  *  This function is experimental and currently only used in polymeres calculation (PEGs or polymeres inside the vesicle).
34  *
35  *  @param *bond is a pointer to a bond between two vertices in polymere
36  *  @param *poly is a pointer to polymere in which we calculate te energy of the bond
37  *  @returns TS_SUCCESS on successful calculation
38 */
fedf2b 39 inline ts_bool bond_energy(ts_bond *bond,ts_poly *poly){
304510 40 //TODO: This value to be changed and implemented in data structure:
M 41     ts_double d_relaxed=1.0;
42     bond->energy=poly->k*pow(bond->bond_length-d_relaxed,2);
fedf2b 43     return TS_SUCCESS;
M 44 };
45
a00f10 46 /** @brief Calculation of energy of the vertex
SP 47  *  
48  *  Main function that calculates energy of the vertex \f$i\f$. Nearest neighbors (NN) must be ordered in counterclockwise direction for this function to work.
49  *  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$.
50  *  
854cb6 51 \begin{tikzpicture}{
a00f10 52 \coordinate[label=below:$i$] (i) at (2,0);
SP 53 \coordinate[label=left:$j_m$] (jm) at (0,3.7);
54 \coordinate[label=above:$j$] (j) at (2.5,6.4);
55 \coordinate[label=right:$j_p$] (jp) at (4,2.7);
d7639a 56
a00f10 57 \draw (i) -- (jm) -- (j) -- (jp) -- (i) -- (j);
SP 58
59 \begin{scope}
60 \path[clip] (jm)--(i)--(j);
61 \draw (jm) circle (0.8);
62 \node[right] at (jm) {$\varphi_m$};
63 \end{scope}
64
65 \begin{scope}
66 \path[clip] (jp)--(i)--(j);
67 \draw (jp) circle (0.8);
68 \node[left] at (jp) {$\varphi_p$};
69 \end{scope}
70
71 %%vertices
72 \draw [fill=gray] (i) circle (0.1);
73 \draw [fill=white] (j) circle (0.1);
74 \draw [fill=white] (jp) circle (0.1);
75 \draw [fill=white] (jm) circle (0.1);
76 %\node[draw,circle,fill=white] at (i) {};
854cb6 77 \end{tikzpicture}
a00f10 78
SP 79  * 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).
80  *
81  * From the curvature the enery is calculated by equation \f$E=\frac{1}{2}\mathbf{h}\cdot\mathbf{h}\f$.
82  * @param *vtx is a pointer to vertex at which we want to calculate the energy
83  * @returns TS_SUCCESS on successful calculation.
84 */
d7639a 85 inline ts_bool energy_vertex(ts_vertex *vtx){
SP 86     ts_uint jj;
87     ts_uint jjp,jjm;
88     ts_vertex *j,*jp, *jm;
89     ts_triangle *jt;
a63f17 90     ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0;
d7639a 91     ts_double x1,x2,x3,ctp,ctm,tot,xlen;
SP 92     ts_double h,ht;
8f6a69 93     for(jj=1; jj<=vtx->neigh_no;jj++){
d7639a 94         jjp=jj+1;
8f6a69 95         if(jjp>vtx->neigh_no) jjp=1;
d7639a 96         jjm=jj-1;
8f6a69 97         if(jjm<1) jjm=vtx->neigh_no;
SP 98         j=vtx->neigh[jj-1];
99         jp=vtx->neigh[jjp-1];
100         jm=vtx->neigh[jjm-1];
101         jt=vtx->tristar[jj-1];
f74313 102         x1=vtx_distance_sq(vtx,jp); //shouldn't be zero!
SP 103         x2=vtx_distance_sq(j,jp); // shouldn't be zero!
8f6a69 104         x3=(j->x-jp->x)*(vtx->x-jp->x)+
SP 105            (j->y-jp->y)*(vtx->y-jp->y)+
106            (j->z-jp->z)*(vtx->z-jp->z);
d7639a 107         
SP 108 #ifdef TS_DOUBLE_DOUBLE
109         ctp=x3/sqrt(x1*x2-x3*x3);
110 #endif
111 #ifdef TS_DOUBLE_FLOAT
112         ctp=x3/sqrtf(x1*x2-x3*x3);
113 #endif
114 #ifdef TS_DOUBLE_LONGDOUBLE
115         ctp=x3/sqrtl(x1*x2-x3*x3);
116 #endif
f74313 117         x1=vtx_distance_sq(vtx,jm);
SP 118         x2=vtx_distance_sq(j,jm);
8f6a69 119         x3=(j->x-jm->x)*(vtx->x-jm->x)+
SP 120            (j->y-jm->y)*(vtx->y-jm->y)+
121            (j->z-jm->z)*(vtx->z-jm->z);
d7639a 122 #ifdef TS_DOUBLE_DOUBLE
SP 123         ctm=x3/sqrt(x1*x2-x3*x3);
124 #endif
125 #ifdef TS_DOUBLE_FLOAT
126         ctm=x3/sqrtf(x1*x2-x3*x3);
127 #endif
128 #ifdef TS_DOUBLE_LONGDOUBLE
129         ctm=x3/sqrtl(x1*x2-x3*x3);
130 #endif
131         tot=ctp+ctm;
132         tot=0.5*tot;
a63f17 133
f74313 134         xlen=vtx_distance_sq(j,vtx);
a63f17 135 /*
d7639a 136 #ifdef  TS_DOUBLE_DOUBLE 
8f6a69 137         vtx->bond[jj-1]->bond_length=sqrt(xlen); 
d7639a 138 #endif
SP 139 #ifdef  TS_DOUBLE_FLOAT
8f6a69 140         vtx->bond[jj-1]->bond_length=sqrtf(xlen); 
d7639a 141 #endif
SP 142 #ifdef  TS_DOUBLE_LONGDOUBLE 
8f6a69 143         vtx->bond[jj-1]->bond_length=sqrtl(xlen); 
d7639a 144 #endif
SP 145
8f6a69 146         vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length;
a63f17 147 */
d7639a 148         s+=tot*xlen;
8f6a69 149         xh+=tot*(j->x - vtx->x);
SP 150         yh+=tot*(j->y - vtx->y);
151         zh+=tot*(j->z - vtx->z);
41a035 152         txn+=jt->xnorm;
SP 153         tyn+=jt->ynorm;
154         tzn+=jt->znorm;
d7639a 155     }
SP 156     
157     h=xh*xh+yh*yh+zh*zh;
158     ht=txn*xh+tyn*yh + tzn*zh;
159     s=s/4.0; 
160 #ifdef TS_DOUBLE_DOUBLE
161     if(ht>=0.0) {
8f6a69 162         vtx->curvature=sqrt(h);
d7639a 163     } else {
8f6a69 164         vtx->curvature=-sqrt(h);
d7639a 165     }
SP 166 #endif
167 #ifdef TS_DOUBLE_FLOAT
168     if(ht>=0.0) {
8f6a69 169         vtx->curvature=sqrtf(h);
d7639a 170     } else {
8f6a69 171         vtx->curvature=-sqrtf(h);
d7639a 172     }
SP 173 #endif
174 #ifdef TS_DOUBLE_LONGDOUBLE
175     if(ht>=0.0) {
8f6a69 176         vtx->curvature=sqrtl(h);
d7639a 177     } else {
8f6a69 178         vtx->curvature=-sqrtl(h);
d7639a 179     }
SP 180 #endif
dac2e5 181 // c is forced curvature energy for each vertex. Should be set to zero for
8f6a69 182 // normal circumstances.
SP 183     vtx->energy=0.5*s*(vtx->curvature/s-vtx->c)*(vtx->curvature/s-vtx->c);
d7639a 184
SP 185     return TS_SUCCESS;
186 }