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