Trisurf Monte Carlo simulator
Samo Penic
2014-03-05 719c9febac2eaff9483fda487b57684afbb59bb2
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
fedf2b 16 ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double *rn){
aec47d 17     ts_uint i;
SP 18     ts_double dist;
19     ts_bool retval; 
20     ts_uint cellidx; 
21     ts_double delta_energy,oenergy;
ed31fe 22     ts_double costheta,sintheta,phi,r;
1ad6d1 23     //This will hold all the information of vtx and its neighbours
dcd350 24     ts_vertex backupvtx[20];
SP 25     memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
a63f17 26
SP 27     //Some stupid tests for debugging cell occupation!
28 /*         cellidx=vertex_self_avoidance(vesicle, vtx);
29     if(vesicle->clist->cell[cellidx]==vtx->cell){
30         fprintf(stderr,"Idx match!\n");
31     } else {
32         fprintf(stderr,"***** Idx don't match!\n");
33         fatal("ENding.",1);
34     }
35 */
36
352fad 37         //temporarly moving the vertex
672ae4 38 //    vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0);
SP 39 //        vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
40 //        vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);
41
ed31fe 42     //random move in a sphere with radius stepsize:
M 43     r=vesicle->stepsize*rn[0];
44     phi=rn[1]*2*M_PI;
45     costheta=2*rn[2]-1;
46     sintheta=sqrt(1-pow(costheta,2));
672ae4 47     vtx->x=vtx->x+r*sintheta*cos(phi);
SP 48     vtx->y=vtx->y+r*sintheta*sin(phi);
49     vtx->z=vtx->z+r*costheta;
50
51
a63f17 52         //distance with neighbours check
6fa0c9 53     for(i=0;i<vtx->neigh->n;i++){
SP 54         dist=vtx_distance_sq(vtx,vtx->neigh->vtx[i]);
8f6a69 55         if(dist<1.0 || dist>vesicle->dmax) {
dcd350 56         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
8f6a69 57         return TS_FAIL;
SP 58         }
aec47d 59     }
304510 60
M 61 // Distance with grafted poly-vertex check:    
62     if(vtx->grafted_poly!=NULL){
63         dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
64         if(dist<1.0 || dist>vesicle->dmax) {
65         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
66         return TS_FAIL;
67         }
68     }
69
aec47d 70     //self avoidance check with distant vertices
352fad 71      cellidx=vertex_self_avoidance(vesicle, vtx);
aec47d 72     //check occupation number
a63f17 73      retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
SP 74     
aec47d 75     if(retval==TS_FAIL){
dcd350 76         vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
aec47d 77         return TS_FAIL;
SP 78     } 
1ad6d1 79    
SP 80  
352fad 81     //if all the tests are successful, then energy for vtx and neighbours is calculated
6fa0c9 82     for(i=0;i<vtx->neigh->n;i++){
SP 83     memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh->vtx[i],sizeof(ts_vertex));
1ad6d1 84     }
aec47d 85
a63f17 86
SP 87
aec47d 88     delta_energy=0;
SP 89     //update the normals of triangles that share bead i.
8f6a69 90     for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17 91     oenergy=vtx->energy;
aec47d 92     energy_vertex(vtx);
a63f17 93     delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d 94     //the same is done for neighbouring vertices
6fa0c9 95     for(i=0;i<vtx->neigh->n;i++){
SP 96         oenergy=vtx->neigh->vtx[i]->energy;
97         energy_vertex(vtx->neigh->vtx[i]);
98         delta_energy+=vtx->neigh->vtx[i]->xk*(vtx->neigh->vtx[i]->energy-oenergy);
aec47d 99     }
304510 100 /* No poly-bond energy for now!
fedf2b 101     if(vtx->grafted_poly!=NULL){
M 102         delta_energy+=
103             (pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
104             pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
105     }
304510 106 */
314f2d 107 //   fprintf(stderr, "DE=%f\n",delta_energy);
aec47d 108     //MONTE CARLOOOOOOOO
SP 109     if(delta_energy>=0){
110 #ifdef TS_DOUBLE_DOUBLE
111         if(exp(-delta_energy)< drand48() )
112 #endif
113 #ifdef TS_DOUBLE_FLOAT
114         if(expf(-delta_energy)< (ts_float)drand48())
115 #endif
116 #ifdef TS_DOUBLE_LONGDOUBLE
117         if(expl(-delta_energy)< (ts_ldouble)drand48())
118 #endif
119     {
120     //not accepted, reverting changes
dcd350 121     vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
6fa0c9 122     for(i=0;i<vtx->neigh->n;i++){
SP 123         vtx->neigh->vtx[i]=memcpy((void *)vtx->neigh->vtx[i],(void *)&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1 124     }
SP 125     
aec47d 126     //update the normals of triangles that share bead i.
dcd350 127    for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
1ad6d1 128
aec47d 129     return TS_FAIL; 
SP 130     }
131 }
a63f17 132         
SP 133 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
134     if(vtx->cell!=vesicle->clist->cell[cellidx]){
135         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
136 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
137         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
138         
139     }
140 //    if(oldcellidx);
aec47d 141     //END MONTE CARLOOOOOOO
SP 142     return TS_SUCCESS;
143 }
144
fedf2b 145
M 146 ts_bool single_poly_vertex_move(ts_vesicle *vesicle,ts_poly *poly,ts_vertex *vtx,ts_double *rn){
147     ts_uint i;
148     ts_bool retval; 
149     ts_uint cellidx; 
304510 150 //    ts_double delta_energy;
fedf2b 151     ts_double costheta,sintheta,phi,r;
304510 152     ts_double dist;
fedf2b 153     //This will hold all the information of vtx and its neighbours
M 154     ts_vertex backupvtx;
304510 155 //    ts_bond backupbond[2];
fedf2b 156     memcpy((void *)&backupvtx,(void *)vtx,sizeof(ts_vertex));
M 157
158     //random move in a sphere with radius stepsize:
159     r=vesicle->stepsize*rn[0];
160     phi=rn[1]*2*M_PI;
161     costheta=2*rn[2]-1;
162     sintheta=sqrt(1-pow(costheta,2));
163     vtx->x=vtx->x+r*sintheta*cos(phi);
164     vtx->y=vtx->y+r*sintheta*sin(phi);
165     vtx->z=vtx->z+r*costheta;
166
167
168     //distance with neighbours check
304510 169     for(i=0;i<vtx->neigh_no;i++){
M 170         dist=vtx_distance_sq(vtx,vtx->neigh[i]);
171         if(dist<1.0 || dist>vesicle->dmax) {
172             vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
173             return TS_FAIL;
174         }
175     }
176
177 // Distance with grafted vesicle-vertex check:    
178     if(vtx==poly->vlist->vtx[0]){
179         dist=vtx_distance_sq(vtx,poly->grafted_vtx);
180         if(dist<1.0 || dist>vesicle->dmax) {
181         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
182         return TS_FAIL;
183         }
184     }
185
fedf2b 186
M 187     //self avoidance check with distant vertices
188     cellidx=vertex_self_avoidance(vesicle, vtx);
189     //check occupation number
190     retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
191     
192     if(retval==TS_FAIL){
193         vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
194         return TS_FAIL;
195     } 
196
197
198     //if all the tests are successful, then energy for vtx and neighbours is calculated
304510 199 /* Energy ignored for now!
fedf2b 200     delta_energy=0;
M 201     for(i=0;i<vtx->bond_no;i++){
202         memcpy((void *)&backupbond[i],(void *)vtx->bond[i],sizeof(ts_bond));
203
204         vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
205         bond_energy(vtx->bond[i],poly);
206         delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
207     }
208
209     if(vtx==poly->vlist->vtx[0]){
210         delta_energy+=
211             (pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
212             pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
213         
214     }
215
216
217     if(delta_energy>=0){
218 #ifdef TS_DOUBLE_DOUBLE
219         if(exp(-delta_energy)< drand48() )
220 #endif
221 #ifdef TS_DOUBLE_FLOAT
222         if(expf(-delta_energy)< (ts_float)drand48())
223 #endif
224 #ifdef TS_DOUBLE_LONGDOUBLE
225         if(expl(-delta_energy)< (ts_ldouble)drand48())
226 #endif
227         {
228     //not accepted, reverting changes
229     vtx=memcpy((void *)vtx,(void *)&backupvtx,sizeof(ts_vertex));
230     for(i=0;i<vtx->bond_no;i++){
231     vtx->bond[i]=memcpy((void *)vtx->bond[i],(void *)&backupbond[i],sizeof(ts_bond));
232     }
233
234     return TS_FAIL; 
235     }
236     }
304510 237 */
fedf2b 238         
M 239 //    oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
240     if(vtx->cell!=vesicle->clist->cell[cellidx]){
241         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
242 //        if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
243         if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);    
244     }
245 //    if(oldcellidx);
246     //END MONTE CARLOOOOOOO
247     return TS_SUCCESS;
248 }