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