Trisurf Monte Carlo simulator
Samo Penic
2019-03-08 7522b886602380764ce76f1d2244fad11bd1054f
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; 
7522b8 23     ts_double delta_energy, oenergy, darea=0.0, dstretchenergy=0.0;
ed31fe 24     ts_double costheta,sintheta,phi,r;
1ad6d1 25     //This will hold all the information of vtx and its neighbours
2afc2f 26     ts_vertex backupvtx[20]; // *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL;
b5cd8c 27     memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
672ae4 28
b5cd8c 29     //random move in a sphere with radius stepsize:
SP 30     r=vesicle->stepsize*rn[0];
31     phi=rn[1]*2*M_PI;
32     costheta=2*rn[2]-1;
33     sintheta=sqrt(1-pow(costheta,2));
34     vtx->x=vtx->x+r*sintheta*cos(phi);
35     vtx->y=vtx->y+r*sintheta*sin(phi);
36     vtx->z=vtx->z+r*costheta;
672ae4 37
fe24d2 38
51b4f0 39 /* Entry point for plugin vm_hard_constraint() function */
77a2c7 40     vesicle->plist->pointer=vesicle->plist->chain->vm_hard_constraint;
SP 41     while(vesicle->plist->pointer!=NULL){
42         retval = vesicle->plist->pointer->plugin->function->vm_hard_constraint(vesicle,vtx, &backupvtx[0]);
51b4f0 43         if(retval==TS_FAIL){
SP 44             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
45             return TS_FAIL;
46         }
77a2c7 47         vesicle->plist->pointer=vesicle->plist->pointer->next;
51b4f0 48     }
36dba8 49 /* End of vm_hard_constraint() */
51b4f0 50
36dba8 51 /* Backuping the neighbours */ 
1ad6d1 52     for(i=0;i<vtx->neigh_no;i++){
dcd350 53     memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
1ad6d1 54     }
aec47d 55
36dba8 56 /* Entry point for plugin vm_energy_before_prepare() */
SP 57
58     vesicle->plist->pointer=vesicle->plist->chain->vm_energy_before_prepare;
59     while(vesicle->plist->pointer!=NULL){
60         vesicle->plist->pointer->plugin->function->vm_energy_before_prepare(vesicle, vtx);
61         vesicle->plist->pointer=vesicle->plist->pointer->next;
c0ae90 62     }
SP 63
64     if(vesicle->tape->constareaswitch==2){
65         for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
66     
67     }
7ec6fb 68     //stretching energy 1 of 3
SP 69     if(vesicle->tape->stretchswitch==1){
699ac4 70         for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
7ec6fb 71     }
aec47d 72     delta_energy=0;
0dd5ba 73
SP 74
fe5069 75 //    vesicle_volume(vesicle);
SP 76 //    fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
43c042 77
aec47d 78     //update the normals of triangles that share bead i.
8f6a69 79     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17 80     oenergy=vtx->energy;
aec47d 81     energy_vertex(vtx);
a63f17 82     delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d 83     //the same is done for neighbouring vertices
8f6a69 84     for(i=0;i<vtx->neigh_no;i++){
SP 85         oenergy=vtx->neigh[i]->energy;
86         energy_vertex(vtx->neigh[i]);
87         delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d 88     }
414b8a 89
77a2c7 90
SP 91 /* Entry point for plugin vm_energy_after_execute() */
92
93     vesicle->plist->pointer=vesicle->plist->chain->vm_energy_after_execute;
94     while(vesicle->plist->pointer!=NULL){
95         delta_energy+=vesicle->plist->pointer->plugin->function->vm_energy_after_execute(vesicle, vtx);
96         vesicle->plist->pointer=vesicle->plist->pointer->next;
97     }
98
99
43c042 100
c0ae90 101     if(vesicle->tape->constareaswitch==2){
SP 102         /* check whether the darea is gt epsarea */
103         for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
104         if(fabs(vesicle->area+darea-A0)>epsarea){
105             //restore old state.
106              vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
107                 for(i=0;i<vtx->neigh_no;i++){
108                     vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
109                 }
110                     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
111                     //fprintf(stderr,"fajlam!\n");
112                     return TS_FAIL;
113         }
114
115
116     }
1121fa 117
250de4 118 /* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */
SP 119     delta_energy+=direct_force_energy(vesicle,vtx,backupvtx);
7ec6fb 120
SP 121     //stretching energy 2 of 3
122     if(vesicle->tape->stretchswitch==1){
123         for(i=0;i<vtx->tristar_no;i++){ 
124             stretchenergy(vesicle, vtx->tristar[i]);
699ac4 125             dstretchenergy+=vtx->tristar[i]->energy;
7ec6fb 126             }
SP 127     }
128
129     delta_energy+=dstretchenergy;    
130         
304510 131 /* No poly-bond energy for now!
fedf2b 132     if(vtx->grafted_poly!=NULL){
M 133         delta_energy+=
134             (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
135             pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
136     }
304510 137 */
0dd5ba 138
SP 139 // plane confinement energy due to compressing force
140     if(vesicle->tape->plane_confinement_switch){
141         if(vesicle->confinement_plane.force_switch){
142             //substract old energy
143             if(abs(vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_max)>1e-10) {
144                 delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-backupvtx[0].z,2);
145                 delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-vtx->z,2);
146             }
147             if(abs(-vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_min)>1e-10) {
148                 delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-backupvtx[0].z,2);
149                 delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-vtx->z,2);
150             }
151         }
152     }
153
77a2c7 154
SP 155 /* Entry point for plugin vm_before_montecarlo_constraint() function */
156     vesicle->plist->pointer=vesicle->plist->chain->vm_before_montecarlo_constraint;
157     while(vesicle->plist->pointer!=NULL){
158         retval = vesicle->plist->pointer->plugin->function->vm_before_montecarlo_constraint(vesicle,vtx, &backupvtx[0]);
159         if(retval==TS_FAIL){
160             vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
161             for(i=0;i<vtx->neigh_no;i++){
162                 vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
163                 }
164             for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); 
165             return TS_FAIL;
166         }
167         vesicle->plist->pointer=vesicle->plist->pointer->next;
168     }
169 /* End of vm_before_montecarlo_constraint() */
170
171
172
173
174
314f2d 175 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 176     //MONTE CARLOOOOOOOO
SP 177     if(delta_energy>=0){
178 #ifdef TS_DOUBLE_DOUBLE
3de289 179         if(exp(-delta_energy)< drand48())
aec47d 180 #endif
SP 181 #ifdef TS_DOUBLE_FLOAT
182         if(expf(-delta_energy)< (ts_float)drand48())
183 #endif
184 #ifdef TS_DOUBLE_LONGDOUBLE
185         if(expl(-delta_energy)< (ts_ldouble)drand48())
186 #endif
187     {
7522b8 188 /*************************************************** MC step rejected **************************************************************/
dcd350 189     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 190     for(i=0;i<vtx->neigh_no;i++){
a63f17 191         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 192     }
SP 193     
7522b8 194     //update the normals of triangles that share bead i.
SP 195     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
196
7ec6fb 197     //stretching energy 3 of 3
SP 198     if(vesicle->tape->stretchswitch==1){
199         for(i=0;i<vtx->tristar_no;i++){ 
200             stretchenergy(vesicle,vtx->tristar[i]);
201             }
202     }
1ad6d1 203
7522b8 204
SP 205
206 /* Entry point for plugin vm_before_montecarlo_constraint() function */
207     vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_rejected;
208     while(vesicle->plist->pointer!=NULL){
209         vesicle->plist->pointer->plugin->function->vm_new_state_rejected(vesicle,vtx, &backupvtx[0]);
210         vesicle->plist->pointer=vesicle->plist->pointer->next;
211     }
212 /* End of vm_before_montecarlo_constraint() */
213
214
aec47d 215     return TS_FAIL; 
SP 216     }
217 }
7522b8 218 /*************************************************** MC step accepted **************************************************************/
51b4f0 219     cellidx=vertex_self_avoidance(vesicle, vtx);
a63f17 220     if(vtx->cell!=vesicle->clist->cell[cellidx]){
SP 221         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
7522b8 222         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);    
a63f17 223     }
2afc2f 224
c0ae90 225
SP 226     if(vesicle->tape->constareaswitch==2){
227         vesicle->area+=darea;
228     }
2afc2f 229
SP 230
231 /* Entry point for plugin vm_before_montecarlo_constraint() function */
232     vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_accepted;
233     while(vesicle->plist->pointer!=NULL){
234         vesicle->plist->pointer->plugin->function->vm_new_state_accepted(vesicle,vtx, &backupvtx[0]);
235         vesicle->plist->pointer=vesicle->plist->pointer->next;
236     }
237 /* End of vm_before_montecarlo_constraint() */
238
aec47d 239     return TS_SUCCESS;
SP 240 }
241
fedf2b 242
7522b8 243
SP 244
245
fedf2b 246 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
M 247     ts_uint i;
248     ts_bool retval; 
249     ts_uint cellidx; 
304510 250 //    ts_double delta_energy;
fedf2b 251     ts_double costheta,sintheta,phi,r;
304510 252     ts_double dist;
fedf2b 253     //This will hold all the information of vtx and its neighbours
M 254     ts_vertex backupvtx;
304510 255 //    ts_bond backupbond[2];
fedf2b 256     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 257
258     //random move in a sphere with radius stepsize:
259     r=vesicle->stepsize*rn[0];
260     phi=rn[1]*2*M_PI;
261     costheta=2*rn[2]-1;
262     sintheta=sqrt(1-pow(costheta,2));
263     vtx->x=vtx->x+r*sintheta*cos(phi);
264     vtx->y=vtx->y+r*sintheta*sin(phi);
265     vtx->z=vtx->z+r*costheta;
266
267
268     //distance with neighbours check
304510 269     for(i=0;i<vtx->neigh_no;i++){
M 270         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
271         if(dist<1.0 || dist>vesicle->dmax) {
272             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
273             return TS_FAIL;
274         }
275     }
276
277 // Distance with grafted vesicle-vertex check:    
278     if(vtx==poly->vlist->vtx[0]){
279         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
280         if(dist<1.0 || dist>vesicle->dmax) {
281         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
282         return TS_FAIL;
283         }
284     }
285
fedf2b 286
M 287     //self avoidance check with distant vertices
288     cellidx=vertex_self_avoidance(vesicle, vtx);
289     //check occupation number
290     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
291     
292     if(retval==TS_FAIL){
293         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
294         return TS_FAIL;
295     } 
296
297
298     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 299 /* Energy ignored for now!
fedf2b 300     delta_energy=0;
M 301     for(i=0;i<vtx->bond_no;i++){
302         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
303
304         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
305         bond_energy(vtx->bond[i],poly);
306         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
307     }
308
309     if(vtx==poly->vlist->vtx[0]){
310         delta_energy+=
311             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
312             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
313         
314     }
315
316
317     if(delta_energy>=0){
318 #ifdef TS_DOUBLE_DOUBLE
319         if(exp(-delta_energy)< drand48() )
320 #endif
321 #ifdef TS_DOUBLE_FLOAT
322         if(expf(-delta_energy)< (ts_float)drand48())
323 #endif
324 #ifdef TS_DOUBLE_LONGDOUBLE
325         if(expl(-delta_energy)< (ts_ldouble)drand48())
326 #endif
327         {
328     //not accepted, reverting changes
329     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
330     for(i=0;i<vtx->bond_no;i++){
331     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
332     }
333
334     return TS_FAIL; 
335     }
336     }
304510 337 */
fedf2b 338         
M 339 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
340     if(vtx->cell!=vesicle->clist->cell[cellidx]){
341         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
342 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
343         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
344     }
345 //    if(oldcellidx);
346     //END MONTE CARLOOOOOOO
347     return TS_SUCCESS;
348 }
58230a 349
M 350
351
352
353 ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
354     ts_uint i;
355     ts_bool retval; 
356     ts_uint cellidx; 
b30f45 357     ts_double delta_energy;
58230a 358     ts_double costheta,sintheta,phi,r;
M 359     ts_double dist[2];
360     //This will hold all the information of vtx and its neighbours
b30f45 361     ts_vertex backupvtx,backupneigh[2];
58230a 362     ts_bond backupbond[2];
b30f45 363
M 364     //backup vertex:        
58230a 365     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 366
367     //random move in a sphere with radius stepsize:
368     r=vesicle->stepsize*rn[0];
369     phi=rn[1]*2*M_PI;
370     costheta=2*rn[2]-1;
371     sintheta=sqrt(1-pow(costheta,2));
372     vtx->x=vtx->x+r*sintheta*cos(phi);
373     vtx->y=vtx->y+r*sintheta*sin(phi);
374     vtx->z=vtx->z+r*costheta;
375
376
377     //distance with neighbours check
378     for(i=0;i<vtx->bond_no;i++){
379         dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
380         if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
381             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
382             return TS_FAIL;
383         }
384     }
385
fe24d2 386 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 387 // Nucleus penetration check:
388     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
389         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
390         return TS_FAIL;
391     }
392
58230a 393
M 394     //self avoidance check with distant vertices
395     cellidx=vertex_self_avoidance(vesicle, vtx);
396     //check occupation number
397     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
398     if(retval==TS_FAIL){
399         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
400         return TS_FAIL;
401     } 
402
403     //backup bonds
404     for(i=0;i<vtx->bond_no;i++){
405         memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
406         vtx->bond[i]->bond_length=sqrt(dist[i]);
407         bond_vector(vtx->bond[i]);
b30f45 408     }
M 409
410     //backup neighboring vertices:
411     for(i=0;i<vtx->neigh_no;i++){
412         memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a 413     }
M 414     
415     //if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45 416     delta_energy=0;
M 417     
418     if(vtx->bond_no == 2){
419         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;
420         delta_energy += vtx->energy - backupvtx.energy;
58230a 421     }
M 422
b30f45 423     for(i=0;i<vtx->neigh_no;i++){
M 424         if(vtx->neigh[i]->bond_no == 2){
425             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;
426             delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
427         }
58230a 428     }
M 429
b30f45 430     // poly->k is filament persistence length (in units l_min)
M 431     delta_energy *= poly->k;
58230a 432
M 433     if(delta_energy>=0){
434 #ifdef TS_DOUBLE_DOUBLE
435         if(exp(-delta_energy)< drand48() )
436 #endif
437 #ifdef TS_DOUBLE_FLOAT
438         if(expf(-delta_energy)< (ts_float)drand48())
439 #endif
440 #ifdef TS_DOUBLE_LONGDOUBLE
441         if(expl(-delta_energy)< (ts_ldouble)drand48())
442 #endif
443         {
444     //not accepted, reverting changes
445     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
b30f45 446     for(i=0;i<vtx->neigh_no;i++){
M 447         memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
448     }
58230a 449     for(i=0;i<vtx->bond_no;i++){
b30f45 450         vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
58230a 451     }
M 452
453     return TS_FAIL; 
454     }
455     }
456     
b30f45 457     
58230a 458 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M 459     if(vtx->cell!=vesicle->clist->cell[cellidx]){
460         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
461 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
462         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
463     }
464 //    if(oldcellidx);
465     //END MONTE CARLOOOOOOO
466     return TS_SUCCESS;
467 }