| 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; | |
| 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 | } |
| SP | 61 | //self avoidance check with distant vertices |
| 352fad | 62 | cellidx=vertex_self_avoidance(vesicle, vtx); |
| aec47d | 63 | //check occupation number |
| a63f17 | 64 | retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx); |
| SP | 65 | |
| aec47d | 66 | if(retval==TS_FAIL){ |
| dcd350 | 67 | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| aec47d | 68 | return TS_FAIL; |
| SP | 69 | } |
| 1ad6d1 | 70 | |
| SP | 71 | |
| 352fad | 72 | //if all the tests are successful, then energy for vtx and neighbours is calculated |
| 1ad6d1 | 73 | for(i=0;i<vtx->neigh_no;i++){ |
| dcd350 | 74 | memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex)); |
| 1ad6d1 | 75 | } |
| aec47d | 76 | |
| a63f17 | 77 | |
| SP | 78 | |
| aec47d | 79 | delta_energy=0; |
| SP | 80 | //update the normals of triangles that share bead i. |
| 8f6a69 | 81 | for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); |
| a63f17 | 82 | oenergy=vtx->energy; |
| aec47d | 83 | energy_vertex(vtx); |
| a63f17 | 84 | delta_energy=vtx->xk*(vtx->energy - oenergy); |
| aec47d | 85 | //the same is done for neighbouring vertices |
| 8f6a69 | 86 | for(i=0;i<vtx->neigh_no;i++){ |
| SP | 87 | oenergy=vtx->neigh[i]->energy; |
| 88 | energy_vertex(vtx->neigh[i]); | |
| 89 | delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy); | |
| aec47d | 90 | } |
| 314f2d | 91 | // fprintf(stderr, "DE=%f\n",delta_energy); |
| aec47d | 92 | //MONTE CARLOOOOOOOO |
| SP | 93 | if(delta_energy>=0){ |
| 94 | #ifdef TS_DOUBLE_DOUBLE | |
| 95 | if(exp(-delta_energy)< drand48() ) | |
| 96 | #endif | |
| 97 | #ifdef TS_DOUBLE_FLOAT | |
| 98 | if(expf(-delta_energy)< (ts_float)drand48()) | |
| 99 | #endif | |
| 100 | #ifdef TS_DOUBLE_LONGDOUBLE | |
| 101 | if(expl(-delta_energy)< (ts_ldouble)drand48()) | |
| 102 | #endif | |
| 103 | { | |
| 104 | //not accepted, reverting changes | |
| dcd350 | 105 | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| 1ad6d1 | 106 | for(i=0;i<vtx->neigh_no;i++){ |
| a63f17 | 107 | vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex)); |
| 1ad6d1 | 108 | } |
| SP | 109 | |
| aec47d | 110 | //update the normals of triangles that share bead i. |
| dcd350 | 111 | for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); |
| 1ad6d1 | 112 | |
| aec47d | 113 | return TS_FAIL; |
| SP | 114 | } |
| 115 | } | |
| a63f17 | 116 | |
| SP | 117 | // oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]); |
| 118 | if(vtx->cell!=vesicle->clist->cell[cellidx]){ | |
| 119 | retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx); | |
| 120 | // if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx); | |
| 121 | if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx); | |
| 122 | ||
| 123 | } | |
| 124 | // if(oldcellidx); | |
| aec47d | 125 | //END MONTE CARLOOOOOOO |
| SP | 126 | return TS_SUCCESS; |
| 127 | } | |
| 128 | ||
