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