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