| commit | author | age | ||
| 7f6076 | 1 | /* vim: set ts=4 sts=4 sw=4 noet : */ |
| d7639a | 2 | #include<stdlib.h> |
| SP | 3 | #include<stdio.h> |
| aec47d | 4 | #include<math.h> |
| SP | 5 | //#include "io.h" |
| 6 | #include "general.h" | |
| 7 | #include "timestep.h" | |
| 8 | #include "vertexmove.h" | |
| 30ee9c | 9 | #include "bondflip.h" |
| d7a113 | 10 | #include "frame.h" |
| SP | 11 | #include "io.h" |
| 37d14a | 12 | #include "stats.h" |
| dc77e8 | 13 | #include "sh.h" |
| 459ff9 | 14 | #include "shcomplex.h" |
| dc77e8 | 15 | #include "vesicle.h" |
| 5a3862 | 16 | #include<gsl/gsl_complex.h> |
| M | 17 | #include<gsl/gsl_complex_math.h> |
| 267db5 | 18 | #include<string.h> |
| ac9826 | 19 | #include <sys/stat.h> |
| 300618 | 20 | #include "initial_distribution.h" |
| fedf2b | 21 | |
| 626811 | 22 | ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){ |
| d5d78a | 23 | ts_uint i, j,k; //,l,m; |
| SP | 24 | ts_double kc1=0,kc2=0,kc3=0,kc4=0; |
| c0ae90 | 25 | ts_double l1,l2,l3,vmsr,bfsr, vmsrt, bfsrt; |
| 37d14a | 26 | ts_ulong epochtime; |
| 0dd5ba | 27 | ts_double max_z,min_z; |
| d5d78a | 28 | FILE *fd3=NULL; |
| 267db5 | 29 | char filename[10000]; |
| ac9826 | 30 | //struct stat st; |
| SP | 31 | strcpy(filename,command_line_args.path); |
| 32 | strcat(filename,"statistics.csv"); | |
| 33 | //int result = stat(filename, &st); | |
| 34 | FILE *fd; | |
| 35 | if(start_iteration==0) | |
| 36 | fd=fopen(filename,"w"); | |
| 37 | else | |
| 38 | fd=fopen(filename,"a"); | |
| 37d14a | 39 | if(fd==NULL){ |
| SP | 40 | fatal("Cannot open statistics.csv file for writing",1); |
| 41 | } | |
| ac9826 | 42 | if(start_iteration==0) |
| SP | 43 | fprintf(fd, "Epoch OuterLoop VertexMoveSucessRate BondFlipSuccessRate Volume Area lamdba1 lambda2 lambda3 Kc(2-9) Kc(6-9) Kc(2-end) Kc(3-6)\n"); |
| 5a3862 | 44 | |
| d5d78a | 45 | /* if(vesicle->sphHarmonics!=NULL){ |
| 267db5 | 46 | strcpy(filename,command_line_args.path); |
| SP | 47 | strcat(filename,"ulm2.csv"); |
| ac9826 | 48 | // int result = stat(filename, &st); |
| SP | 49 | if(start_iteration==0) |
| 267db5 | 50 | fd2=fopen(filename,"w"); |
| ac9826 | 51 | else |
| SP | 52 | fd2=fopen(filename,"a"); |
| 819a09 | 53 | if(fd2==NULL){ |
| S | 54 | fatal("Cannot open ulm2.csv file for writing",1); |
| d5d78a | 55 | } |
| SP | 56 | |
| 819a09 | 57 | if(start_iteration==0) //file does not exist |
| S | 58 | fprintf(fd2, "Timestep u_00^2 u_10^2 u_11^2 u_20^2 ...\n"); |
| 5a3862 | 59 | } |
| d5d78a | 60 | */ |
| 5a3862 | 61 | |
| c60a49 | 62 | /* RANDOM SEED SET BY CURRENT TIME */ |
| M | 63 | epochtime=get_epoch(); |
| 64 | srand48(epochtime); | |
| d7a113 | 65 | centermass(vesicle); |
| SP | 66 | cell_occupation(vesicle); |
| fe5069 | 67 | vesicle_volume(vesicle); //needed for constant volume at this moment |
| 88bdd7 | 68 | vesicle_area(vesicle); //needed for constant area at this moment |
| 49981c | 69 | if(V0<0.000001) |
| SP | 70 | V0=vesicle->volume; |
| 71 | ts_fprintf(stdout,"Setting volume V0=%.17f\n",V0); | |
| 72 | if(A0<0.000001) | |
| 73 | A0=vesicle->area; | |
| 819a09 | 74 | ts_fprintf(stdout,"Setting area A0=%.17f\n",A0); |
| a54977 | 75 | epsvol=4.0*sqrt(2.0*M_PI)/pow(3.0,3.0/4.0)*V0/pow(vesicle->tlist->n,3.0/2.0); |
| d5d78a | 76 | // printf("epsvol=%e\n",epsvol); |
| 88bdd7 | 77 | epsarea=A0/(ts_double)vesicle->tlist->n; |
| 300618 | 78 | |
| SP | 79 | int delta_number_c0=ceil((float)vesicle->tape->number_of_vertices_with_c0/(float)(iterations-1)); |
| 80 | ts_fprintf(stdout, "Adding additional %d inclusions every %d iterations (%d times, total=%d).\n",delta_number_c0, mcsweeps, (vesicle->tape->number_of_vertices_with_c0)/delta_number_c0,(vesicle->tape->number_of_vertices_with_c0/delta_number_c0)*delta_number_c0); | |
| 88bdd7 | 81 | |
| 626811 | 82 | 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 | 83 | for(i=start_iteration;i<inititer+iterations;i++){ |
| 37d14a | 84 | vmsr=0.0; |
| SP | 85 | bfsr=0.0; |
| 0dd5ba | 86 | |
| SP | 87 | //plane confinement |
| 88 | if(vesicle->tape->plane_confinement_switch){ | |
| 89 | min_z=1e10; | |
| 90 | max_z=-1e10; | |
| 91 | for(k=0;k<vesicle->vlist->n;k++){ | |
| 92 | if(vesicle->vlist->vtx[k]->z > max_z) max_z=vesicle->vlist->vtx[k]->z; | |
| 93 | if(vesicle->vlist->vtx[k]->z < min_z) min_z=vesicle->vlist->vtx[k]->z; | |
| 94 | } | |
| 95 | vesicle->confinement_plane.force_switch=0; | |
| 96 | if(max_z>=vesicle->tape->plane_d/2.0){ | |
| 97 | ts_fprintf(stdout, "Max vertex out of bounds (z>=%e). Plane set to max_z = %e.\n",vesicle->tape->plane_d/2.0,max_z); | |
| 98 | vesicle->confinement_plane.z_max = max_z; | |
| 99 | vesicle->confinement_plane.force_switch=1; | |
| 100 | } else { | |
| 101 | vesicle->confinement_plane.z_max=vesicle->tape->plane_d/2.0; | |
| 102 | } | |
| 103 | if(min_z<=-vesicle->tape->plane_d/2.0){ | |
| 104 | ts_fprintf(stdout, "Min vertex out of bounds (z<=%e). Plane set to min_z = %e.\n",-vesicle->tape->plane_d/2.0,min_z); | |
| 105 | vesicle->confinement_plane.z_min = min_z; | |
| 106 | vesicle->confinement_plane.force_switch=1; | |
| 107 | } else { | |
| 108 | vesicle->confinement_plane.z_min=-vesicle->tape->plane_d/2.0; | |
| 109 | } | |
| 110 | ts_fprintf(stdout,"Vesicle confinement by plane set to (zmin, zmax)=(%e,%e).\n",vesicle->confinement_plane.z_min,vesicle->confinement_plane.z_max); | |
| 111 | if(vesicle->confinement_plane.force_switch) ts_fprintf(stdout,"Squeezing with force %e.\n",vesicle->tape->plane_F); | |
| 112 | } | |
| 113 | ||
| 114 | //end plane confinement | |
| 115 | ||
| 3de289 | 116 | /* vesicle_volume(vesicle); |
| SP | 117 | fprintf(stderr,"Volume before TS=%1.16e\n", vesicle->volume); */ |
| 300618 | 118 | |
| SP | 119 | if(i*delta_number_c0<vesicle->tape->number_of_vertices_with_c0){ |
| 120 | fprintf(stdout, "Adding %d inclusions.\n", delta_number_c0); | |
| 121 | add_vertices_with_c0(vesicle, delta_number_c0, vesicle->tape->c0, vesicle->tape->w); | |
| 122 | } | |
| d7a113 | 123 | for(j=0;j<mcsweeps;j++){ |
| 37d14a | 124 | single_timestep(vesicle, &vmsrt, &bfsrt); |
| SP | 125 | vmsr+=vmsrt; |
| 126 | bfsr+=bfsrt; | |
| d7a113 | 127 | } |
| 3de289 | 128 | /* |
| SP | 129 | vesicle_volume(vesicle); |
| 130 | fprintf(stderr,"Volume after TS=%1.16e\n", vesicle->volume); */ | |
| 37d14a | 131 | vmsr/=(ts_double)mcsweeps; |
| SP | 132 | bfsr/=(ts_double)mcsweeps; |
| d7a113 | 133 | centermass(vesicle); |
| SP | 134 | cell_occupation(vesicle); |
| 2a1e3d | 135 | /* BINARY DUMPS ARE OBSOLETE. SHOULD WORK AS OF MAR 2020, BUT NO LONGER MAINTAINED */ |
| SP | 136 | // dump_state(vesicle,i); |
| d5d78a | 137 | vesicle_volume(vesicle); //calculates just volume. |
| SP | 138 | vesicle_area(vesicle); //calculates area. |
| 139 | if(vesicle->tape->constvolswitch==0){ | |
| 140 | V0=vesicle->volume; | |
| 141 | } | |
| 142 | if(vesicle->tape->constareaswitch==0){ | |
| 143 | A0=vesicle->area; | |
| 144 | } | |
| 58230a | 145 | if(i>=inititer){ |
| 0a2c81 | 146 | write_vertex_xml_file(vesicle,i-inititer,NULL); |
| 267db5 | 147 | write_master_xml_file(command_line_args.output_fullfilename); |
| 37d14a | 148 | epochtime=get_epoch(); |
| SP | 149 | gyration_eigen(vesicle, &l1, &l2, &l3); |
| d5d78a | 150 | //r0=getR0(vesicle); |
| dc77e8 | 151 | |
| c0ae90 | 152 | fprintf(fd, "%lu %u %e %e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e %1.16e\n",epochtime,i,vmsr,bfsr,vesicle->volume, vesicle->area,l1,l2,l3,kc1, kc2, kc3,kc4); |
| 5a3862 | 153 | |
| 632960 | 154 | fflush(fd); |
| 144784 | 155 | // sprintf(filename,"timestep-%05d.pov",i-inititer); |
| fe24d2 | 156 | // write_pov_file(vesicle,filename); |
| 49981c | 157 | } //end if(inititer....) |
| SP | 158 | fd3=fopen(".status","w"); //write status file when everything is written to disk. |
| 159 | if(fd3==NULL){ | |
| 160 | fatal("Cannot open .status file for writing",1); | |
| d7a113 | 161 | } |
| 49981c | 162 | fprintf(fd3,"%d",i); |
| SP | 163 | fclose(fd3); |
| 164 | ts_fprintf(stdout,"Done %d out of %d iterations (x %d MC sweeps).\n",i+1,inititer+iterations,mcsweeps); | |
| d7a113 | 165 | } |
| 37d14a | 166 | fclose(fd); |
| d5d78a | 167 | // if(fd2!=NULL) fclose(fd2); |
| d7a113 | 168 | return TS_SUCCESS; |
| SP | 169 | } |
| d7639a | 170 | |
| 37d14a | 171 | ts_bool single_timestep(ts_vesicle *vesicle,ts_double *vmsr, ts_double *bfsr){ |
| 3de289 | 172 | // vesicle_volume(vesicle); |
| SP | 173 | // fprintf(stderr,"Volume before TS=%1.16e\n", vesicle->volume); |
| d7639a | 174 | ts_bool retval; |
| SP | 175 | ts_double rnvec[3]; |
| fe5069 | 176 | ts_uint i,j, b; |
| 37d14a | 177 | ts_uint vmsrcnt=0; |
| aec47d | 178 | for(i=0;i<vesicle->vlist->n;i++){ |
| d7639a | 179 | rnvec[0]=drand48(); |
| SP | 180 | rnvec[1]=drand48(); |
| 181 | rnvec[2]=drand48(); | |
| aec47d | 182 | retval=single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec); |
| 37d14a | 183 | if(retval==TS_SUCCESS) vmsrcnt++; |
| d7639a | 184 | } |
| SP | 185 | |
| 37d14a | 186 | ts_int bfsrcnt=0; |
| fedf2b | 187 | for(i=0;i<3*vesicle->vlist->n;i++){ |
| fe5069 | 188 | b=rand() % vesicle->blist->n; |
| d7639a | 189 | //find a bond and return a pointer to a bond... |
| SP | 190 | //call single_bondflip_timestep... |
| fe5069 | 191 | retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec); |
| 3de289 | 192 | // b++; retval=TS_FAIL; |
| 37d14a | 193 | if(retval==TS_SUCCESS) bfsrcnt++; |
| fedf2b | 194 | } |
| M | 195 | |
| 196 | for(i=0;i<vesicle->poly_list->n;i++){ | |
| 58230a | 197 | for(j=0;j<vesicle->poly_list->poly[i]->vlist->n;j++){ |
| M | 198 | rnvec[0]=drand48(); |
| 199 | rnvec[1]=drand48(); | |
| 200 | rnvec[2]=drand48(); | |
| 201 | retval=single_poly_vertex_move(vesicle,vesicle->poly_list->poly[i],vesicle->poly_list->poly[i]->vlist->vtx[j],rnvec); | |
| 202 | } | |
| fedf2b | 203 | } |
| M | 204 | |
| 58230a | 205 | |
| M | 206 | for(i=0;i<vesicle->filament_list->n;i++){ |
| 207 | for(j=0;j<vesicle->filament_list->poly[i]->vlist->n;j++){ | |
| 208 | rnvec[0]=drand48(); | |
| 209 | rnvec[1]=drand48(); | |
| 210 | rnvec[2]=drand48(); | |
| 211 | retval=single_filament_vertex_move(vesicle,vesicle->filament_list->poly[i],vesicle->filament_list->poly[i]->vlist->vtx[j],rnvec); | |
| 212 | } | |
| fedf2b | 213 | } |
| M | 214 | |
| 58230a | 215 | |
| fedf2b | 216 | // printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,3*vesicle->blist->n,(double)cnt/(double)vesicle->blist->n/3.0); |
| 37d14a | 217 | *vmsr=(ts_double)vmsrcnt/(ts_double)vesicle->vlist->n; |
| SP | 218 | *bfsr=(ts_double)bfsrcnt/(ts_double)vesicle->vlist->n/3.0; |
| 3de289 | 219 | // vesicle_volume(vesicle); |
| SP | 220 | // fprintf(stderr,"Volume after TS=%1.16e\n", vesicle->volume); |
| d7639a | 221 | return TS_SUCCESS; |
| SP | 222 | } |
| 223 | ||
| 224 | ||
| 225 | ||
