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