Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2019-03-09 9a1f16ae3affc4db83f2eb2623418ac5cff2af03
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
aec47d 2 #include<stdlib.h>
SP 3 #include<math.h>
4 #include "general.h"
5 #include "vertex.h"
6 #include "bond.h"
7 #include "triangle.h"
8 #include "vesicle.h"
9 #include "energy.h"
10 #include "timestep.h"
11 #include "cell.h"
9166cb 12 #include "io.h"
aec47d 13 #include<stdio.h>
SP 14 #include "vertexmove.h"
1ad6d1 15 #include <string.h>
43c042 16 #include "constvol.h"
51b4f0 17 #include "plugins.h"
aec47d 18
fedf2b 19 ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
aec47d 20     ts_uint i;
SP 21     ts_bool retval; 
22     ts_uint cellidx; 
2c4278 23     ts_double delta_energy, oenergy;
ed31fe 24     ts_double costheta,sintheta,phi,r;
9a1f16 25     ts_vertex backupvtx[20];
b5cd8c 26     memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
672ae4 27
b5cd8c 28     //random move in a sphere with radius stepsize:
SP 29     r=vesicle->stepsize*rn[0];
30     phi=rn[1]*2*M_PI;
31     costheta=2*rn[2]-1;
32     sintheta=sqrt(1-pow(costheta,2));
33     vtx->x=vtx->x+r*sintheta*cos(phi);
34     vtx->y=vtx->y+r*sintheta*sin(phi);
35     vtx->z=vtx->z+r*costheta;
672ae4 36
fe24d2 37
51b4f0 38 /* Entry point for plugin vm_hard_constraint() function */
77a2c7 39     vesicle->plist->pointer=vesicle->plist->chain->vm_hard_constraint;
SP 40     while(vesicle->plist->pointer!=NULL){
41         retval = vesicle->plist->pointer->plugin->function->vm_hard_constraint(vesicle,vtx, &backupvtx[0]);
51b4f0 42         if(retval==TS_FAIL){
SP 43             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
44             return TS_FAIL;
45         }
77a2c7 46         vesicle->plist->pointer=vesicle->plist->pointer->next;
51b4f0 47     }
36dba8 48 /* End of vm_hard_constraint() */
51b4f0 49
36dba8 50 /* Backuping the neighbours */ 
1ad6d1 51     for(i=0;i<vtx->neigh_no;i++){
dcd350 52     memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
1ad6d1 53     }
aec47d 54
36dba8 55 /* Entry point for plugin vm_energy_before_prepare() */
SP 56     vesicle->plist->pointer=vesicle->plist->chain->vm_energy_before_prepare;
57     while(vesicle->plist->pointer!=NULL){
58         vesicle->plist->pointer->plugin->function->vm_energy_before_prepare(vesicle, vtx);
59         vesicle->plist->pointer=vesicle->plist->pointer->next;
c0ae90 60     }
5e08f2 61 /* End of vm_energy_before_prepare() */
c0ae90 62
9a1f16 63     //update the normals of triangles that share bead i.
SP 64     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17 65     oenergy=vtx->energy;
9a1f16 66     energy_vertex(vtx);
SP 67     delta_energy=vtx->xk*(vtx->energy - oenergy);
68     //the same is done for neighbouring vertices
69     for(i=0;i<vtx->neigh_no;i++){
70         oenergy=vtx->neigh[i]->energy;
71         energy_vertex(vtx->neigh[i]);
72         delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
73     }
414b8a 74
77a2c7 75 /* Entry point for plugin vm_energy_after_execute() */
SP 76     vesicle->plist->pointer=vesicle->plist->chain->vm_energy_after_execute;
77     while(vesicle->plist->pointer!=NULL){
85898e 78         delta_energy+=vesicle->plist->pointer->plugin->function->vm_energy_after_execute(vesicle, vtx, backupvtx);
77a2c7 79         vesicle->plist->pointer=vesicle->plist->pointer->next;
SP 80     }
2c4278 81 /* End of vm_energy_after_execute() */
77a2c7 82
SP 83
304510 84 /* No poly-bond energy for now!
fedf2b 85     if(vtx->grafted_poly!=NULL){
M 86         delta_energy+=
87             (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
88             pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
89     }
304510 90 */
0dd5ba 91
77a2c7 92
SP 93 /* Entry point for plugin vm_before_montecarlo_constraint() function */
94     vesicle->plist->pointer=vesicle->plist->chain->vm_before_montecarlo_constraint;
95     while(vesicle->plist->pointer!=NULL){
96         retval = vesicle->plist->pointer->plugin->function->vm_before_montecarlo_constraint(vesicle,vtx, &backupvtx[0]);
97         if(retval==TS_FAIL){
98             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
99             for(i=0;i<vtx->neigh_no;i++){
100                 vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
101                 }
102             for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
103             return TS_FAIL;
104         }
105         vesicle->plist->pointer=vesicle->plist->pointer->next;
106     }
107 /* End of vm_before_montecarlo_constraint() */
108
109
110
111
112
314f2d 113 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 114     //MONTE CARLOOOOOOOO
SP 115     if(delta_energy>=0){
116 #ifdef TS_DOUBLE_DOUBLE
3de289 117         if(exp(-delta_energy)< drand48())
aec47d 118 #endif
SP 119 #ifdef TS_DOUBLE_FLOAT
120         if(expf(-delta_energy)< (ts_float)drand48())
121 #endif
122 #ifdef TS_DOUBLE_LONGDOUBLE
123         if(expl(-delta_energy)< (ts_ldouble)drand48())
124 #endif
125     {
7522b8 126 /*************************************************** MC step rejected **************************************************************/
dcd350 127     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 128     for(i=0;i<vtx->neigh_no;i++){
a63f17 129         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 130     }
SP 131     
7522b8 132     //update the normals of triangles that share bead i.
SP 133     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
134
135
136 /* Entry point for plugin vm_before_montecarlo_constraint() function */
137     vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_rejected;
138     while(vesicle->plist->pointer!=NULL){
139         vesicle->plist->pointer->plugin->function->vm_new_state_rejected(vesicle,vtx, &backupvtx[0]);
140         vesicle->plist->pointer=vesicle->plist->pointer->next;
141     }
142 /* End of vm_before_montecarlo_constraint() */
143
144
aec47d 145     return TS_FAIL; 
SP 146     }
147 }
7522b8 148 /*************************************************** MC step accepted **************************************************************/
51b4f0 149     cellidx=vertex_self_avoidance(vesicle, vtx);
a63f17 150     if(vtx->cell!=vesicle->clist->cell[cellidx]){
SP 151         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
7522b8 152         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);    
a63f17 153     }
2afc2f 154
SP 155 /* Entry point for plugin vm_before_montecarlo_constraint() function */
156     vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_accepted;
157     while(vesicle->plist->pointer!=NULL){
158         vesicle->plist->pointer->plugin->function->vm_new_state_accepted(vesicle,vtx, &backupvtx[0]);
159         vesicle->plist->pointer=vesicle->plist->pointer->next;
160     }
161 /* End of vm_before_montecarlo_constraint() */
162
aec47d 163     return TS_SUCCESS;
SP 164 }
165
fedf2b 166
7522b8 167
SP 168
169
fedf2b 170 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
M 171     ts_uint i;
172     ts_bool retval; 
173     ts_uint cellidx; 
304510 174 //    ts_double delta_energy;
fedf2b 175     ts_double costheta,sintheta,phi,r;
304510 176     ts_double dist;
fedf2b 177     //This will hold all the information of vtx and its neighbours
M 178     ts_vertex backupvtx;
304510 179 //    ts_bond backupbond[2];
fedf2b 180     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 181
182     //random move in a sphere with radius stepsize:
183     r=vesicle->stepsize*rn[0];
184     phi=rn[1]*2*M_PI;
185     costheta=2*rn[2]-1;
186     sintheta=sqrt(1-pow(costheta,2));
187     vtx->x=vtx->x+r*sintheta*cos(phi);
188     vtx->y=vtx->y+r*sintheta*sin(phi);
189     vtx->z=vtx->z+r*costheta;
190
191
192     //distance with neighbours check
304510 193     for(i=0;i<vtx->neigh_no;i++){
M 194         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
195         if(dist<1.0 || dist>vesicle->dmax) {
196             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
197             return TS_FAIL;
198         }
199     }
200
201 // Distance with grafted vesicle-vertex check:    
202     if(vtx==poly->vlist->vtx[0]){
203         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
204         if(dist<1.0 || dist>vesicle->dmax) {
205         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
206         return TS_FAIL;
207         }
208     }
209
fedf2b 210
M 211     //self avoidance check with distant vertices
212     cellidx=vertex_self_avoidance(vesicle, vtx);
213     //check occupation number
214     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
215     
216     if(retval==TS_FAIL){
217         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
218         return TS_FAIL;
219     } 
220
221
222     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 223 /* Energy ignored for now!
fedf2b 224     delta_energy=0;
M 225     for(i=0;i<vtx->bond_no;i++){
226         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
227
228         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
229         bond_energy(vtx->bond[i],poly);
230         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
231     }
232
233     if(vtx==poly->vlist->vtx[0]){
234         delta_energy+=
235             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
236             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
237         
238     }
239
240
241     if(delta_energy>=0){
242 #ifdef TS_DOUBLE_DOUBLE
243         if(exp(-delta_energy)< drand48() )
244 #endif
245 #ifdef TS_DOUBLE_FLOAT
246         if(expf(-delta_energy)< (ts_float)drand48())
247 #endif
248 #ifdef TS_DOUBLE_LONGDOUBLE
249         if(expl(-delta_energy)< (ts_ldouble)drand48())
250 #endif
251         {
252     //not accepted, reverting changes
253     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
254     for(i=0;i<vtx->bond_no;i++){
255     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
256     }
257
258     return TS_FAIL; 
259     }
260     }
304510 261 */
fedf2b 262         
M 263 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
264     if(vtx->cell!=vesicle->clist->cell[cellidx]){
265         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
266 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
267         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
268     }
269 //    if(oldcellidx);
270     //END MONTE CARLOOOOOOO
271     return TS_SUCCESS;
272 }
58230a 273
M 274
275
276
277 ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
278     ts_uint i;
279     ts_bool retval; 
280     ts_uint cellidx; 
b30f45 281     ts_double delta_energy;
58230a 282     ts_double costheta,sintheta,phi,r;
M 283     ts_double dist[2];
284     //This will hold all the information of vtx and its neighbours
b30f45 285     ts_vertex backupvtx,backupneigh[2];
58230a 286     ts_bond backupbond[2];
b30f45 287
M 288     //backup vertex:        
58230a 289     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 290
291     //random move in a sphere with radius stepsize:
292     r=vesicle->stepsize*rn[0];
293     phi=rn[1]*2*M_PI;
294     costheta=2*rn[2]-1;
295     sintheta=sqrt(1-pow(costheta,2));
296     vtx->x=vtx->x+r*sintheta*cos(phi);
297     vtx->y=vtx->y+r*sintheta*sin(phi);
298     vtx->z=vtx->z+r*costheta;
299
300
301     //distance with neighbours check
302     for(i=0;i<vtx->bond_no;i++){
303         dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
304         if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
305             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
306             return TS_FAIL;
307         }
308     }
309
fe24d2 310 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 311 // Nucleus penetration check:
312     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
313         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
314         return TS_FAIL;
315     }
316
58230a 317
M 318     //self avoidance check with distant vertices
319     cellidx=vertex_self_avoidance(vesicle, vtx);
320     //check occupation number
321     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
322     if(retval==TS_FAIL){
323         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
324         return TS_FAIL;
325     } 
326
327     //backup bonds
328     for(i=0;i<vtx->bond_no;i++){
329         memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
330         vtx->bond[i]->bond_length=sqrt(dist[i]);
331         bond_vector(vtx->bond[i]);
b30f45 332     }
M 333
334     //backup neighboring vertices:
335     for(i=0;i<vtx->neigh_no;i++){
336         memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a 337     }
M 338     
339     //if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45 340     delta_energy=0;
M 341     
342     if(vtx->bond_no == 2){
343         vtx->energy = -(vtx->bond[0]->x*vtx->bond[1]->x + vtx->bond[0]->y*vtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
344         delta_energy += vtx->energy - backupvtx.energy;
58230a 345     }
M 346
b30f45 347     for(i=0;i<vtx->neigh_no;i++){
M 348         if(vtx->neigh[i]->bond_no == 2){
349             vtx->neigh[i]->energy = -(vtx->neigh[i]->bond[0]->x*vtx->neigh[i]->bond[1]->x + vtx->neigh[i]->bond[0]->y*vtx->neigh[i]->bond[1]->y + vtx->neigh[i]->bond[0]->z*vtx->neigh[i]->bond[1]->z)/vtx->neigh[i]->bond[0]->bond_length/vtx->neigh[i]->bond[1]->bond_length;
350             delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
351         }
58230a 352     }
M 353
b30f45 354     // poly->k is filament persistence length (in units l_min)
M 355     delta_energy *= poly->k;
58230a 356
M 357     if(delta_energy>=0){
358 #ifdef TS_DOUBLE_DOUBLE
359         if(exp(-delta_energy)< drand48() )
360 #endif
361 #ifdef TS_DOUBLE_FLOAT
362         if(expf(-delta_energy)< (ts_float)drand48())
363 #endif
364 #ifdef TS_DOUBLE_LONGDOUBLE
365         if(expl(-delta_energy)< (ts_ldouble)drand48())
366 #endif
367         {
368     //not accepted, reverting changes
369     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
b30f45 370     for(i=0;i<vtx->neigh_no;i++){
M 371         memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
372     }
58230a 373     for(i=0;i<vtx->bond_no;i++){
b30f45 374         vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
58230a 375     }
M 376
377     return TS_FAIL; 
378     }
379     }
380     
b30f45 381     
58230a 382 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M 383     if(vtx->cell!=vesicle->clist->cell[cellidx]){
384         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
385 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
386         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
387     }
388 //    if(oldcellidx);
389     //END MONTE CARLOOOOOOO
390     return TS_SUCCESS;
391 }