Trisurf Monte Carlo simulator
Samo Penic
2014-03-08 1153c110781fc3add2c22ab2d996e05b0e0b10de
commit | author | age
aec47d 1 #include<stdlib.h>
SP 2 #include<math.h>
3 #include "general.h"
4 #include "vertex.h"
5 #include "bond.h"
6 #include "triangle.h"
7 #include "vesicle.h"
8 #include "energy.h"
9 #include "timestep.h"
10 #include "cell.h"
11 //#include "io.h"
12 #include<stdio.h>
13 #include "vertexmove.h"
1ad6d1 14 #include <string.h>
aec47d 15
fedf2b 16 ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
aec47d 17     ts_uint i;
SP 18     ts_double dist;
19     ts_bool retval; 
20     ts_uint cellidx; 
414b8a 21     ts_double delta_energy,oenergy,dvol=0.0;
ed31fe 22     ts_double costheta,sintheta,phi,r;
1ad6d1 23     //This will hold all the information of vtx and its neighbours
dcd350 24     ts_vertex backupvtx[20];
SP 25     memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
a63f17 26
SP 27     //Some stupid tests for debugging cell occupation!
28 /*         cellidx=vertex_self_avoidance(vesicle, vtx);
29     if(vesicle->clist->cell[cellidx]==vtx->cell){
30         fprintf(stderr,"Idx match!\n");
31     } else {
32         fprintf(stderr,"***** Idx don't match!\n");
33         fatal("ENding.",1);
34     }
35 */
36
352fad 37         //temporarly moving the vertex
672ae4 38 //    vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
SP 39 //        vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
40 //        vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
41
ed31fe 42     //random move in a sphere with radius stepsize:
M 43     r=vesicle->stepsize*rn[0];
44     phi=rn[1]*2*M_PI;
45     costheta=2*rn[2]-1;
46     sintheta=sqrt(1-pow(costheta,2));
672ae4 47     vtx->x=vtx->x+r*sintheta*cos(phi);
SP 48     vtx->y=vtx->y+r*sintheta*sin(phi);
49     vtx->z=vtx->z+r*costheta;
50
51
a63f17 52         //distance with neighbours check
8f6a69 53     for(i=0;i<vtx->neigh_no;i++){
352fad 54         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
8f6a69 55         if(dist<1.0 || dist>vesicle->dmax) {
dcd350 56         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
8f6a69 57         return TS_FAIL;
SP 58         }
aec47d 59     }
304510 60
M 61 // Distance with grafted poly-vertex check:    
62     if(vtx->grafted_poly!=NULL){
63         dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
64         if(dist<1.0 || dist>vesicle->dmax) {
65         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
66         return TS_FAIL;
67         }
68     }
69
aec47d 70     //self avoidance check with distant vertices
352fad 71      cellidx=vertex_self_avoidance(vesicle, vtx);
aec47d 72     //check occupation number
a63f17 73      retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
SP 74     
aec47d 75     if(retval==TS_FAIL){
dcd350 76         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
aec47d 77         return TS_FAIL;
SP 78     } 
1ad6d1 79    
SP 80  
352fad 81     //if all the tests are successful, then energy for vtx and neighbours is calculated
1ad6d1 82     for(i=0;i<vtx->neigh_no;i++){
dcd350 83     memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
1ad6d1 84     }
aec47d 85
414b8a 86     if(vesicle->pswitch == 1){
M 87         for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
88     };
a63f17 89
aec47d 90     delta_energy=0;
SP 91     //update the normals of triangles that share bead i.
8f6a69 92     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17 93     oenergy=vtx->energy;
aec47d 94     energy_vertex(vtx);
a63f17 95     delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d 96     //the same is done for neighbouring vertices
8f6a69 97     for(i=0;i<vtx->neigh_no;i++){
SP 98         oenergy=vtx->neigh[i]->energy;
99         energy_vertex(vtx->neigh[i]);
100         delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d 101     }
414b8a 102
M 103     if(vesicle->pswitch == 1){
104         for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
105         delta_energy-=vesicle->pressure*dvol;
106     };
107
304510 108 /* No poly-bond energy for now!
fedf2b 109     if(vtx->grafted_poly!=NULL){
M 110         delta_energy+=
111             (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
112             pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
113     }
304510 114 */
314f2d 115 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 116     //MONTE CARLOOOOOOOO
SP 117     if(delta_energy>=0){
118 #ifdef TS_DOUBLE_DOUBLE
119         if(exp(-delta_energy)< drand48() )
120 #endif
121 #ifdef TS_DOUBLE_FLOAT
122         if(expf(-delta_energy)< (ts_float)drand48())
123 #endif
124 #ifdef TS_DOUBLE_LONGDOUBLE
125         if(expl(-delta_energy)< (ts_ldouble)drand48())
126 #endif
127     {
128     //not accepted, reverting changes
dcd350 129     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 130     for(i=0;i<vtx->neigh_no;i++){
a63f17 131         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 132     }
SP 133     
aec47d 134     //update the normals of triangles that share bead i.
dcd350 135    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1 136
aec47d 137     return TS_FAIL; 
SP 138     }
139 }
a63f17 140         
SP 141 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
142     if(vtx->cell!=vesicle->clist->cell[cellidx]){
143         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
144 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
145         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
146         
147     }
148 //    if(oldcellidx);
aec47d 149     //END MONTE CARLOOOOOOO
SP 150     return TS_SUCCESS;
151 }
152
fedf2b 153
M 154 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
155     ts_uint i;
156     ts_bool retval; 
157     ts_uint cellidx; 
304510 158 //    ts_double delta_energy;
fedf2b 159     ts_double costheta,sintheta,phi,r;
304510 160     ts_double dist;
fedf2b 161     //This will hold all the information of vtx and its neighbours
M 162     ts_vertex backupvtx;
304510 163 //    ts_bond backupbond[2];
fedf2b 164     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 165
166     //random move in a sphere with radius stepsize:
167     r=vesicle->stepsize*rn[0];
168     phi=rn[1]*2*M_PI;
169     costheta=2*rn[2]-1;
170     sintheta=sqrt(1-pow(costheta,2));
171     vtx->x=vtx->x+r*sintheta*cos(phi);
172     vtx->y=vtx->y+r*sintheta*sin(phi);
173     vtx->z=vtx->z+r*costheta;
174
175
176     //distance with neighbours check
304510 177     for(i=0;i<vtx->neigh_no;i++){
M 178         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
179         if(dist<1.0 || dist>vesicle->dmax) {
180             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
181             return TS_FAIL;
182         }
183     }
184
185 // Distance with grafted vesicle-vertex check:    
186     if(vtx==poly->vlist->vtx[0]){
187         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
188         if(dist<1.0 || dist>vesicle->dmax) {
189         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
190         return TS_FAIL;
191         }
192     }
193
fedf2b 194
M 195     //self avoidance check with distant vertices
196     cellidx=vertex_self_avoidance(vesicle, vtx);
197     //check occupation number
198     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
199     
200     if(retval==TS_FAIL){
201         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
202         return TS_FAIL;
203     } 
204
205
206     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 207 /* Energy ignored for now!
fedf2b 208     delta_energy=0;
M 209     for(i=0;i<vtx->bond_no;i++){
210         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
211
212         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
213         bond_energy(vtx->bond[i],poly);
214         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
215     }
216
217     if(vtx==poly->vlist->vtx[0]){
218         delta_energy+=
219             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
220             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
221         
222     }
223
224
225     if(delta_energy>=0){
226 #ifdef TS_DOUBLE_DOUBLE
227         if(exp(-delta_energy)< drand48() )
228 #endif
229 #ifdef TS_DOUBLE_FLOAT
230         if(expf(-delta_energy)< (ts_float)drand48())
231 #endif
232 #ifdef TS_DOUBLE_LONGDOUBLE
233         if(expl(-delta_energy)< (ts_ldouble)drand48())
234 #endif
235         {
236     //not accepted, reverting changes
237     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
238     for(i=0;i<vtx->bond_no;i++){
239     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
240     }
241
242     return TS_FAIL; 
243     }
244     }
304510 245 */
fedf2b 246         
M 247 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
248     if(vtx->cell!=vesicle->clist->cell[cellidx]){
249         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
250 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
251         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
252     }
253 //    if(oldcellidx);
254     //END MONTE CARLOOOOOOO
255     return TS_SUCCESS;
256 }