Added constant area. Unfortunately dump has changed sue to change in datast...

Samo Penic

2014-09-18 c0ae90d68622063a44af443acc01f7e4925fbc56

Added constant area. Unfortunately dump has changed sue to change in datastructure

 src/bondflip.c

@@ -31,7 +31,7 @@
    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;
@@ -170,7 +170,7 @@
  //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);
*/
@@ -192,6 +192,39 @@
        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 */
@@ -305,12 +338,14 @@
     /* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
//            fprintf(stderr,"SUCCESS!!!\n");

 if(vesicle->tape->constvolswitch == 2){
    vesicle->volume+=dvol;
    } else
    if(vesicle->tape->constvolswitch == 1){
    if(vesicle->tape->constvolswitch == 2){
        vesicle->volume+=dvol;
    } 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);

 src/general.h

@@ -261,6 +261,7 @@
    long int R_nucleus;
    long int pswitch;
    long int constvolswitch;
    long int constareaswitch;
    ts_double constvolprecision;
        char *multiprocessing;
       long int brezveze0;
@@ -305,7 +306,7 @@
    ts_int pswitch;
    ts_tape *tape;
    ts_double R_nucleus;

    ts_double area;
} ts_vesicle;


@@ -314,8 +315,9 @@

int quiet;
ts_double V0;
ts_double A0;
ts_double epsvol;

ts_double epsarea;
/* FUNCTIONS */

/** Non-fatal error function handler:

 src/io.c

@@ -1008,6 +1008,7 @@
        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),

 src/tape

@@ -20,10 +20,12 @@
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

 src/timestep.c

@@ -19,7 +19,7 @@
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];
@@ -49,8 +49,11 @@
    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++){
@@ -77,8 +80,8 @@
            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);
@@ -117,7 +120,7 @@

            }

            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);

 src/vertexmove.c

@@ -20,7 +20,7 @@
    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;
@@ -41,7 +41,7 @@
//        vtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0);
//        vtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0);

    //random move in a sphere with radius stepsize:
//random move in a sphere with radius stepsize:
    r=vesicle->stepsize*rn[0];
    phi=rn[1]*2*M_PI;
    costheta=2*rn[2]-1;
@@ -51,12 +51,12 @@
    vtx->z=vtx->z+r*costheta;


        //distance with neighbours check
//distance with neighbours check
    for(i=0;i<vtx->neigh_no;i++){
        dist=vtx_distance_sq(vtx,vtx->neigh[i]);
        if(dist<1.0 || dist>vesicle->dmax) {
        vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
        return TS_FAIL;
            vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex));
            return TS_FAIL;
        }
    }

@@ -87,14 +87,19 @@
    } 
   
 
    //if all the tests are successful, then energy for vtx and neighbours is calculated
//if all the tests are successful, then energy for vtx and neighbours is calculated
    for(i=0;i<vtx->neigh_no;i++){
    memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
    }

    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;
    
@@ -118,6 +123,22 @@
        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 */
@@ -211,6 +232,10 @@
    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);

 src/vesicle.c

@@ -60,3 +60,20 @@
    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;
}

 src/vesicle.h

@@ -5,5 +5,5 @@
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