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