Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2019-03-08 77a2c7de0dd6563a30467fa2e73f02ec9f4c3deb
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; 
77a2c7 23     ts_double delta_energy, oenergy,dvol=0.0, 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
b5cd8c 26     ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL;
SP 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     {
188     //not accepted, reverting changes
fbcbdf 189   //  fprintf(stderr,"MC failed\n");
dcd350 190     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
1ad6d1 191     for(i=0;i<vtx->neigh_no;i++){
a63f17 192         vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 193     }
SP 194     
aec47d 195     //update the normals of triangles that share bead i.
dcd350 196    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
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
fe5069 204 //    fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
43c042 205     if(vesicle->tape->constvolswitch == 1){
958e0e 206         constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
43c042 207     }
fe5069 208 //    fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
dd5aca 209 //    vesicle_volume(vesicle);
SP 210 //    fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
aec47d 211     return TS_FAIL; 
SP 212     }
213 }
2b14da 214     //accepted    
fbcbdf 215  //   fprintf(stderr,"MC accepted\n");
a63f17 216 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
51b4f0 217     cellidx=vertex_self_avoidance(vesicle, vtx);
a63f17 218     if(vtx->cell!=vesicle->clist->cell[cellidx]){
SP 219         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
220 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
221         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
222         
223     }
2b14da 224
1121fa 225     if(vesicle->tape->constvolswitch == 2){
SP 226     vesicle->volume+=dvol;
227     } else
43c042 228     if(vesicle->tape->constvolswitch == 1){
fbcbdf 229         constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
43c042 230     }
c0ae90 231
SP 232     if(vesicle->tape->constareaswitch==2){
233         vesicle->area+=darea;
234     }
a63f17 235 //    if(oldcellidx);
aec47d 236     //END MONTE CARLOOOOOOO
dd5aca 237 //    vesicle_volume(vesicle);
SP 238 //    fprintf(stderr,"Volume after success=%1.16e\n", vesicle->volume);
aec47d 239     return TS_SUCCESS;
SP 240 }
241
fedf2b 242
M 243 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
244     ts_uint i;
245     ts_bool retval; 
246     ts_uint cellidx; 
304510 247 //    ts_double delta_energy;
fedf2b 248     ts_double costheta,sintheta,phi,r;
304510 249     ts_double dist;
fedf2b 250     //This will hold all the information of vtx and its neighbours
M 251     ts_vertex backupvtx;
304510 252 //    ts_bond backupbond[2];
fedf2b 253     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 254
255     //random move in a sphere with radius stepsize:
256     r=vesicle->stepsize*rn[0];
257     phi=rn[1]*2*M_PI;
258     costheta=2*rn[2]-1;
259     sintheta=sqrt(1-pow(costheta,2));
260     vtx->x=vtx->x+r*sintheta*cos(phi);
261     vtx->y=vtx->y+r*sintheta*sin(phi);
262     vtx->z=vtx->z+r*costheta;
263
264
265     //distance with neighbours check
304510 266     for(i=0;i<vtx->neigh_no;i++){
M 267         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
268         if(dist<1.0 || dist>vesicle->dmax) {
269             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
270             return TS_FAIL;
271         }
272     }
273
274 // Distance with grafted vesicle-vertex check:    
275     if(vtx==poly->vlist->vtx[0]){
276         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
277         if(dist<1.0 || dist>vesicle->dmax) {
278         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
279         return TS_FAIL;
280         }
281     }
282
fedf2b 283
M 284     //self avoidance check with distant vertices
285     cellidx=vertex_self_avoidance(vesicle, vtx);
286     //check occupation number
287     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
288     
289     if(retval==TS_FAIL){
290         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
291         return TS_FAIL;
292     } 
293
294
295     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 296 /* Energy ignored for now!
fedf2b 297     delta_energy=0;
M 298     for(i=0;i<vtx->bond_no;i++){
299         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
300
301         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
302         bond_energy(vtx->bond[i],poly);
303         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
304     }
305
306     if(vtx==poly->vlist->vtx[0]){
307         delta_energy+=
308             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
309             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
310         
311     }
312
313
314     if(delta_energy>=0){
315 #ifdef TS_DOUBLE_DOUBLE
316         if(exp(-delta_energy)< drand48() )
317 #endif
318 #ifdef TS_DOUBLE_FLOAT
319         if(expf(-delta_energy)< (ts_float)drand48())
320 #endif
321 #ifdef TS_DOUBLE_LONGDOUBLE
322         if(expl(-delta_energy)< (ts_ldouble)drand48())
323 #endif
324         {
325     //not accepted, reverting changes
326     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
327     for(i=0;i<vtx->bond_no;i++){
328     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
329     }
330
331     return TS_FAIL; 
332     }
333     }
304510 334 */
fedf2b 335         
M 336 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
337     if(vtx->cell!=vesicle->clist->cell[cellidx]){
338         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
339 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
340         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
341     }
342 //    if(oldcellidx);
343     //END MONTE CARLOOOOOOO
344     return TS_SUCCESS;
345 }
58230a 346
M 347
348
349
350 ts_bool single_filament_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
351     ts_uint i;
352     ts_bool retval; 
353     ts_uint cellidx; 
b30f45 354     ts_double delta_energy;
58230a 355     ts_double costheta,sintheta,phi,r;
M 356     ts_double dist[2];
357     //This will hold all the information of vtx and its neighbours
b30f45 358     ts_vertex backupvtx,backupneigh[2];
58230a 359     ts_bond backupbond[2];
b30f45 360
M 361     //backup vertex:        
58230a 362     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 363
364     //random move in a sphere with radius stepsize:
365     r=vesicle->stepsize*rn[0];
366     phi=rn[1]*2*M_PI;
367     costheta=2*rn[2]-1;
368     sintheta=sqrt(1-pow(costheta,2));
369     vtx->x=vtx->x+r*sintheta*cos(phi);
370     vtx->y=vtx->y+r*sintheta*sin(phi);
371     vtx->z=vtx->z+r*costheta;
372
373
374     //distance with neighbours check
375     for(i=0;i<vtx->bond_no;i++){
376         dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
377         if(dist[i]<1.0 || dist[i]>vesicle->dmax) {
378             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
379             return TS_FAIL;
380         }
381     }
382
fe24d2 383 // TODO: Maybe faster if checks only nucleus-neighboring cells
M 384 // Nucleus penetration check:
385     if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
386         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
387         return TS_FAIL;
388     }
389
58230a 390
M 391     //self avoidance check with distant vertices
392     cellidx=vertex_self_avoidance(vesicle, vtx);
393     //check occupation number
394     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
395     if(retval==TS_FAIL){
396         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
397         return TS_FAIL;
398     } 
399
400     //backup bonds
401     for(i=0;i<vtx->bond_no;i++){
402         memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
403         vtx->bond[i]->bond_length=sqrt(dist[i]);
404         bond_vector(vtx->bond[i]);
b30f45 405     }
M 406
407     //backup neighboring vertices:
408     for(i=0;i<vtx->neigh_no;i++){
409         memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a 410     }
M 411     
412     //if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45 413     delta_energy=0;
M 414     
415     if(vtx->bond_no == 2){
416         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;
417         delta_energy += vtx->energy - backupvtx.energy;
58230a 418     }
M 419
b30f45 420     for(i=0;i<vtx->neigh_no;i++){
M 421         if(vtx->neigh[i]->bond_no == 2){
422             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;
423             delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
424         }
58230a 425     }
M 426
b30f45 427     // poly->k is filament persistence length (in units l_min)
M 428     delta_energy *= poly->k;
58230a 429
M 430     if(delta_energy>=0){
431 #ifdef TS_DOUBLE_DOUBLE
432         if(exp(-delta_energy)< drand48() )
433 #endif
434 #ifdef TS_DOUBLE_FLOAT
435         if(expf(-delta_energy)< (ts_float)drand48())
436 #endif
437 #ifdef TS_DOUBLE_LONGDOUBLE
438         if(expl(-delta_energy)< (ts_ldouble)drand48())
439 #endif
440         {
441     //not accepted, reverting changes
442     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
b30f45 443     for(i=0;i<vtx->neigh_no;i++){
M 444         memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
445     }
58230a 446     for(i=0;i<vtx->bond_no;i++){
b30f45 447         vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
58230a 448     }
M 449
450     return TS_FAIL; 
451     }
452     }
453     
b30f45 454     
58230a 455 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M 456     if(vtx->cell!=vesicle->clist->cell[cellidx]){
457         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
458 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
459         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
460     }
461 //    if(oldcellidx);
462     //END MONTE CARLOOOOOOO
463     return TS_SUCCESS;
464 }