Added constant area. Unfortunately dump has changed sue to change in datastructure
| | |
| | | ts_vertex *k=bond->vtx2; |
| | | ts_uint nei,neip,neim; |
| | | ts_uint i,j; |
| | | ts_double oldenergy, delta_energy, dvol=0.0; |
| | | ts_double oldenergy, delta_energy, dvol=0.0, darea=0.0; |
| | | ts_triangle *lm=NULL,*lp=NULL, *lp1=NULL, *lm2=NULL; |
| | | |
| | | ts_vertex *kp,*km; |
| | |
| | | //Neigbours of k, it, km, kp don't change its energy. |
| | | |
| | | if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){dvol = -lm->volume - lp->volume;} |
| | | |
| | | if(vesicle->tape->constareaswitch==2){darea=-lm->area-lp->area;} |
| | | /* vesicle_volume(vesicle); |
| | | fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume); |
| | | */ |
| | |
| | | dvol = dvol + lm->volume + lp->volume; |
| | | if(vesicle->pswitch==1) delta_energy-= vesicle->pressure*dvol; |
| | | } |
| | | if(vesicle->tape->constareaswitch==2){ |
| | | darea=darea+lm->area+lp->area; |
| | | /*check whether the dvol is gt than epsvol */ |
| | | if(fabs(vesicle->area+darea-A0)>epsarea){ |
| | | //restore old state. |
| | | /* restoration procedure copied from few lines below */ |
| | | for(i=0;i<4;i++){ |
| | | // fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no ); |
| | | free(orig_vtx[i]->neigh); |
| | | free(orig_vtx[i]->tristar); |
| | | free(orig_vtx[i]->bond); |
| | | free(orig_tria[i]->neigh); |
| | | memcpy((void *)orig_vtx[i],(void *)bck_vtx[i],sizeof(ts_vertex)); |
| | | memcpy((void *)orig_tria[i],(void *)bck_tria[i],sizeof(ts_triangle)); |
| | | // fprintf(stderr,"Restored vtx neigh[%d] with neighbours %d\n",i, orig_vtx[i]->neigh_no ); |
| | | /* level 2 pointers are redirected*/ |
| | | } |
| | | memcpy(bond,bck_bond,sizeof(ts_bond)); |
| | | for(i=0;i<4;i++){ |
| | | free(bck_vtx[i]); |
| | | free(bck_tria[i]); |
| | | /* fprintf(stderr,"Restoring vtx neigh[%d] with neighbours %d =",i, orig_vtx[i]->neigh_no ); |
| | | for(j=0;j<orig_vtx[i]->neigh_no;j++) fprintf(stderr," %d", orig_vtx[i]->neigh[j]->idx); |
| | | fprintf(stderr,"\n"); */ |
| | | } |
| | | free(bck_bond); |
| | | return TS_FAIL; |
| | | |
| | | } |
| | | } |
| | | |
| | | |
| | | |
| | | |
| | | if(vesicle->tape->constvolswitch == 2){ |
| | | /*check whether the dvol is gt than epsvol */ |
| | |
| | | |
| | | if(vesicle->tape->constvolswitch == 2){ |
| | | vesicle->volume+=dvol; |
| | | } else |
| | | if(vesicle->tape->constvolswitch == 1){ |
| | | } else if(vesicle->tape->constvolswitch == 1){ |
| | | constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup); |
| | | } |
| | | if(vesicle->tape->constareaswitch==2){ |
| | | vesicle->area+=darea; |
| | | } |
| | | // delete all backups |
| | | for(i=0;i<4;i++){ |
| | | free(bck_vtx[i]->neigh); |
| | |
| | | long int R_nucleus; |
| | | long int pswitch; |
| | | long int constvolswitch; |
| | | long int constareaswitch; |
| | | ts_double constvolprecision; |
| | | char *multiprocessing; |
| | | long int brezveze0; |
| | |
| | | ts_int pswitch; |
| | | ts_tape *tape; |
| | | ts_double R_nucleus; |
| | | |
| | | ts_double area; |
| | | } ts_vesicle; |
| | | |
| | | |
| | |
| | | |
| | | int quiet; |
| | | ts_double V0; |
| | | ts_double A0; |
| | | ts_double epsvol; |
| | | |
| | | ts_double epsarea; |
| | | /* FUNCTIONS */ |
| | | |
| | | /** Non-fatal error function handler: |
| | |
| | | CFG_SIMPLE_FLOAT("xk0",&tape->xk0), |
| | | CFG_SIMPLE_INT("pswitch",&tape->pswitch), |
| | | CFG_SIMPLE_INT("constvolswitch",&tape->constvolswitch), |
| | | CFG_SIMPLE_INT("constareaswitch",&tape->constareaswitch), |
| | | CFG_SIMPLE_FLOAT("constvolprecision",&tape->constvolprecision), |
| | | CFG_SIMPLE_FLOAT("pressure",&tape->pressure), |
| | | CFG_SIMPLE_FLOAT("k_spring",&tape->kspring), |
| | |
| | | constvolswitch=2 |
| | | constvolprecision=1e-14 |
| | | |
| | | #Constant area constraint (0 disable constant area, 2 enable constant area with epsarea) |
| | | constareaswitch=2 |
| | | |
| | | ####### Polymer (brush) definitions ########### |
| | | # npoly is a number of polymers attached to npoly distinct vertices on vesicle |
| | | npoly=5 |
| | | npoly=0 |
| | | # nmono is a number of monomers in each polymer |
| | | nmono=10 |
| | | # Spring constant between monomers of the polymer |
| | |
| | | ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){ |
| | | ts_uint i, j,k,l,m; |
| | | ts_double r0,kc1,kc2,kc3,kc4; |
| | | ts_double l1,l2,l3,volume=0.0,area=0.0,vmsr,bfsr, vmsrt, bfsrt; |
| | | ts_double l1,l2,l3,vmsr,bfsr, vmsrt, bfsrt; |
| | | ts_ulong epochtime; |
| | | FILE *fd1,*fd2=NULL; |
| | | char filename[10000]; |
| | |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | vesicle_volume(vesicle); //needed for constant volume at this moment |
| | | vesicle_area(vesicle); //needed for constant area at this moment |
| | | V0=vesicle->volume; |
| | | A0=vesicle->area; |
| | | epsvol=4.0*sqrt(2.0*M_PI)/pow(3.0,3.0/4.0)*V0/pow(vesicle->tlist->n,3.0/2.0); |
| | | epsarea=A0/(ts_double)vesicle->tlist->n; |
| | | // fprintf(stderr, "DVol=%1.16f (%1.16f), V0=%1.16f\n", epsvol,0.003e-2*V0,V0); |
| | | if(start_iteration<inititer) ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps); |
| | | for(i=start_iteration;i<inititer+iterations;i++){ |
| | |
| | | write_master_xml_file(command_line_args.output_fullfilename); |
| | | epochtime=get_epoch(); |
| | | gyration_eigen(vesicle, &l1, &l2, &l3); |
| | | vesicle_volume(vesicle); //calculates just volume. Area is not added to ts_vesicle yet! |
| | | get_area_volume(vesicle, &area,&volume); //that's why I must recalculate area (and volume for no particular reason). |
| | | vesicle_volume(vesicle); //calculates just volume. |
| | | vesicle_area(vesicle); //calculates area. |
| | | r0=getR0(vesicle); |
| | | if(vesicle->sphHarmonics!=NULL){ |
| | | preparationSh(vesicle,r0); |
| | |
| | | |
| | | } |
| | | |
| | | 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,volume, area,l1,l2,l3,kc1, kc2, kc3,kc4); |
| | | 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); |
| | | |
| | | fflush(fd); |
| | | // sprintf(filename,"timestep-%05d.pov",i-inititer); |
| | |
| | | ts_double dist; |
| | | ts_bool retval; |
| | | ts_uint cellidx; |
| | | ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0; |
| | | ts_double delta_energy, delta_energy_cv,oenergy,dvol=0.0, darea=0.0; |
| | | ts_double costheta,sintheta,phi,r; |
| | | //This will hold all the information of vtx and its neighbours |
| | | ts_vertex backupvtx[20], *constvol_vtx_moved=NULL, *constvol_vtx_backup=NULL; |
| | |
| | | |
| | | if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){ |
| | | for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume; |
| | | }; |
| | | } |
| | | |
| | | if(vesicle->tape->constareaswitch==2){ |
| | | for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area; |
| | | |
| | | } |
| | | |
| | | delta_energy=0; |
| | | |
| | |
| | | if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol; |
| | | }; |
| | | |
| | | if(vesicle->tape->constareaswitch==2){ |
| | | /* check whether the darea is gt epsarea */ |
| | | for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area; |
| | | if(fabs(vesicle->area+darea-A0)>epsarea){ |
| | | //restore old state. |
| | | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| | | for(i=0;i<vtx->neigh_no;i++){ |
| | | vtx->neigh[i]=memcpy((void *)vtx->neigh[i],(void *)&backupvtx[i+1],sizeof(ts_vertex)); |
| | | } |
| | | for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); |
| | | //fprintf(stderr,"fajlam!\n"); |
| | | return TS_FAIL; |
| | | } |
| | | |
| | | |
| | | } |
| | | |
| | | if(vesicle->tape->constvolswitch==2){ |
| | | /*check whether the dvol is gt than epsvol */ |
| | |
| | | if(vesicle->tape->constvolswitch == 1){ |
| | | constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup); |
| | | } |
| | | |
| | | if(vesicle->tape->constareaswitch==2){ |
| | | vesicle->area+=darea; |
| | | } |
| | | // if(oldcellidx); |
| | | //END MONTE CARLOOOOOOO |
| | | // vesicle_volume(vesicle); |
| | |
| | | vesicle->volume=volume; |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | | /* @brief Function makes a sum of partial areas of each triangle. |
| | | * |
| | | * |
| | | * |
| | | */ |
| | | ts_bool vesicle_area(ts_vesicle *vesicle){ |
| | | ts_double area; |
| | | ts_uint i; |
| | | ts_triangle **tria=vesicle->tlist->tria; |
| | | area=0; |
| | | for(i=0;i<vesicle->tlist->n;i++){ |
| | | area=area+tria[i]->area; |
| | | } |
| | | vesicle->area=area; |
| | | return TS_SUCCESS; |
| | | } |
| | |
| | | ts_bool vesicle_translate(ts_vesicle *vesicle,ts_double x, ts_double y, ts_double z); |
| | | ts_bool vesicle_free(ts_vesicle *vesicle); |
| | | ts_bool vesicle_volume(ts_vesicle *vesicle); |
| | | |
| | | ts_bool vesicle_area(ts_vesicle *vesicle); |
| | | #endif |