Trisurf Monte Carlo simulator
mihaf
2014-03-25 fe24d29b3c8684f08dccd01f5785aa48b7137322
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include<stdio.h>
aec47d 3 #include<math.h>
SP 4 //#include "io.h"
5 #include "general.h"
6 #include "timestep.h"
7 #include "vertexmove.h"
30ee9c 8 #include "bondflip.h"
d7a113 9 #include "frame.h"
SP 10 #include "io.h"
fedf2b 11
626811 12 ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
d7a113 13     ts_uint i, j;
fe24d2 14 //     char filename[255];
d7a113 15
SP 16     centermass(vesicle);
17     cell_occupation(vesicle);
626811 18     if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
SP 19     for(i=start_iteration;i<inititer+iterations;i++){
d7a113 20         for(j=0;j<mcsweeps;j++){
SP 21             single_timestep(vesicle);
22         }
23         centermass(vesicle);
24         cell_occupation(vesicle);
f8e6ba 25         ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps);
1ab449 26             dump_state(vesicle,i);
58230a 27         if(i>=inititer){
d7a113 28             write_vertex_xml_file(vesicle,i-inititer);
fe24d2 29     write_master_xml_file("test.pvd");
M 30
144784 31         //    sprintf(filename,"timestep-%05d.pov",i-inititer);
fe24d2 32         //    write_pov_file(vesicle,filename);
d7a113 33         }
SP 34     }
35     return TS_SUCCESS;
36 }
d7639a 37
SP 38 ts_bool single_timestep(ts_vesicle *vesicle){
39     ts_bool retval;
40     ts_double rnvec[3];
fedf2b 41     ts_uint i,j,b;
aec47d 42     for(i=0;i<vesicle->vlist->n;i++){
d7639a 43         rnvec[0]=drand48();
SP 44         rnvec[1]=drand48();
45         rnvec[2]=drand48();
aec47d 46         retval=single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec);
d7639a 47     }
SP 48
b7c9b3 49 //    ts_int cnt=0;
fedf2b 50     for(i=0;i<3*vesicle->vlist->n;i++){
dcf17d 51 //why is rnvec needed in bondflip?
SP 52 /*        rnvec[0]=drand48();
53         rnvec[1]=drand48();
54         rnvec[2]=drand48();
55 */ 
142a67 56     b=rand() % vesicle->blist->n;
d7639a 57         //find a bond and return a pointer to a bond...
SP 58         //call single_bondflip_timestep...
142a67 59         retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
b7c9b3 60 //    if(retval==TS_SUCCESS) cnt++;        
fedf2b 61     }
M 62
63     for(i=0;i<vesicle->poly_list->n;i++){
58230a 64         for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){
M 65             rnvec[0]=drand48();
66             rnvec[1]=drand48();
67             rnvec[2]=drand48();
68             retval=single_poly_vertex_move(vesicle,vesicle->poly_list->poly[i],vesicle->poly_list->poly[i]->vlist->vtx[j],rnvec);    
69         }
fedf2b 70     }
M 71
58230a 72
M 73     for(i=0;i<vesicle->filament_list->n;i++){
74         for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){
75             rnvec[0]=drand48();
76             rnvec[1]=drand48();
77             rnvec[2]=drand48();
78             retval=single_filament_vertex_move(vesicle,vesicle->filament_list->poly[i],vesicle->filament_list->poly[i]->vlist->vtx[j],rnvec);    
79         }
fedf2b 80     }
M 81  
58230a 82
fedf2b 83 //    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,3*vesicle->blist->n,(double)cnt/(double)vesicle->blist->n/3.0);
41a035 84     if(retval);
d7639a 85     return TS_SUCCESS;
SP 86 }
87
88
89