| 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" | |
| 14 | ||
| 15 | ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double | |
| 16 | *rn){ | |
| 17 | ts_uint i; | |
| 18 | ts_double dist; | |
| 19 | ts_vertex *tvtx=(ts_vertex *)malloc(sizeof(ts_vertex)); | |
| 20 | tvtx->data=init_vertex_data(); | |
| 21 | ts_bool retval; | |
| 22 | ts_uint cellidx; | |
| 23 | ts_double xold,yold,zold; | |
| 24 | ts_double delta_energy,oenergy; | |
| 25 | ts_vertex *ovtx; | |
| 26 | ||
| 27 | //randomly we move the temporary vertex | |
| 28 | tvtx->data->x=vtx->data->x+vesicle->stepsize*(2.0*rn[0]-1.0); | |
| 29 | tvtx->data->y=vtx->data->y+vesicle->stepsize*(2.0*rn[1]-1.0); | |
| 30 | tvtx->data->z=vtx->data->z+vesicle->stepsize*(2.0*rn[2]-1.0); | |
| 31 | //check we if some length to neighbours are too much | |
| 32 | for(i=0;i<vtx->data->neigh_no;i++){ | |
| 33 | dist=vtx_distance_sq(tvtx,vtx->data->neigh[i]); | |
| 34 | if(dist<1.0 || dist>vesicle->dmax) return TS_FAIL; | |
| 35 | } | |
| 36 | fprintf(stderr,"Was here!\n"); | |
| 37 | //self avoidance check with distant vertices | |
| 38 | cellidx=vertex_self_avoidance(vesicle, tvtx); | |
| 39 | //check occupation number | |
| 40 | retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx); | |
| 41 | if(retval==TS_FAIL){ | |
| 42 | return TS_FAIL; | |
| 43 | } | |
| 44 | ||
| 45 | //if all the tests are successful, then we update the vertex position | |
| 46 | xold=vtx->data->x; | |
| 47 | yold=vtx->data->y; | |
| 48 | zold=vtx->data->z; | |
| 49 | ovtx=malloc(sizeof(ts_vertex)); | |
| 50 | vtx_copy(ovtx,vtx); | |
| 51 | vtx->data->x=tvtx->data->x; | |
| 52 | vtx->data->y=tvtx->data->y; | |
| 53 | vtx->data->z=tvtx->data->z; | |
| 54 | ||
| 55 | delta_energy=0; | |
| 56 | //update the normals of triangles that share bead i. | |
| 57 | for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]); | |
| 58 | //energy and curvature | |
| 59 | energy_vertex(vtx); | |
| 60 | delta_energy=vtx->data->xk*(vtx->data->energy - ovtx->data->energy); | |
| 61 | //the same is done for neighbouring vertices | |
| 62 | for(i=0;i<vtx->data->neigh_no;i++){ | |
| 63 | oenergy=vtx->data->neigh[i]->data->energy; | |
| 64 | energy_vertex(vtx->data->neigh[i]); | |
| 65 | delta_energy+=vtx->data->neigh[i]->data->xk*(vtx->data->neigh[i]->data->energy-oenergy); | |
| 66 | } | |
| 67 | fprintf(stderr, "DE=%f\n",delta_energy); | |
| 68 | //MONTE CARLOOOOOOOO | |
| 69 | if(delta_energy>=0){ | |
| 70 | #ifdef TS_DOUBLE_DOUBLE | |
| 71 | if(exp(-delta_energy)< drand48() ) | |
| 72 | #endif | |
| 73 | #ifdef TS_DOUBLE_FLOAT | |
| 74 | if(expf(-delta_energy)< (ts_float)drand48()) | |
| 75 | #endif | |
| 76 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 77 | if(expl(-delta_energy)< (ts_ldouble)drand48()) | |
| 78 | #endif | |
| 79 | { | |
| 80 | //not accepted, reverting changes | |
| 81 | vtx->data->x=xold; | |
| 82 | vtx->data->y=yold; | |
| 83 | vtx->data->z=zold; | |
| 84 | //update the normals of triangles that share bead i. | |
| 85 | for(i=0;i<vtx->data->tristar_no;i++) triangle_normal_vector(vtx->data->tristar[i]); | |
| 86 | //energy and curvature | |
| 87 | energy_vertex(vtx); | |
| 88 | //the same is done for neighbouring vertices | |
| 89 | for(i=0;i<vtx->data->neigh_no;i++) energy_vertex(vtx->data->neigh[i]); | |
| 90 | free(ovtx->data->bond_length); | |
| 91 | free(ovtx->data->bond_length_dual); | |
| 92 | free(ovtx); | |
| 93 | return TS_FAIL; | |
| 94 | } | |
| 95 | } | |
| 96 | //END MONTE CARLOOOOOOO | |
| 97 | ||
| 98 | //TODO: change cell occupation if necessary! | |
| 99 | ||
| 100 | free(ovtx->data->bond_length); | |
| 101 | free(ovtx->data->bond_length_dual); | |
| 102 | free(ovtx); | |
| 103 | return TS_SUCCESS; | |
| 104 | } | |
| 105 | ||
