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