First variant of constant volume constrant (a new one proposed by Miha afte...

Samo Penic

2014-09-05 1121faa13a2038facad22073f0fc610903d98691

First variant of constant volume constrant (a new one proposed by Miha after reading his article). It seems to work, however there are still some things to be done, such as Miha's derivation of the epsvol (0.003% is used at the moment) and solving the problem of additional global variables.

 src/bondflip.c

@@ -169,7 +169,7 @@
  oldenergy+=it->xk* it->energy;
  //Neigbours of k, it, km, kp don't change its energy.

    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch==1){dvol = -lm->volume - lp->volume;}
    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){dvol = -lm->volume - lp->volume;}
      
/*    vesicle_volume(vesicle);
    fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
@@ -188,11 +188,41 @@
  //Neigbours of k, it, km, kp don't change its energy.

    delta_energy-=oldenergy;
    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch==1){
    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
        dvol = dvol + lm->volume + lp->volume;
        if(vesicle->pswitch==1) delta_energy-= vesicle->pressure*dvol;
    }

    if(vesicle->tape->constvolswitch == 2){
        /*check whether the dvol is gt than epsvol */
        if(fabs(vesicle->volume+dvol-V0)>epsvol){
            //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;

        }

    } else
    if(vesicle->tape->constvolswitch == 1){
        retval=constvolume(vesicle, it, -dvol, &delta_energy_cv, &constvol_vtx_moved,&constvol_vtx_backup);
        if(retval==TS_FAIL){
@@ -275,6 +305,9 @@
     /* IF BONDFLIP ACCEPTED, THEN RETURN SUCCESS! */
//            fprintf(stderr,"SUCCESS!!!\n");

 if(vesicle->tape->constvolswitch == 2){
    vesicle->volume+=dvol;
    } else
    if(vesicle->tape->constvolswitch == 1){
        constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
    }

 src/general.h

@@ -313,7 +313,8 @@
/* GLOBAL VARIABLES */

int quiet;

ts_double V0;
ts_double epsvol;

/* FUNCTIONS */


 src/timestep.c

@@ -45,6 +45,8 @@
    centermass(vesicle);
    cell_occupation(vesicle);
    vesicle_volume(vesicle); //needed for constant volume at this moment
    V0=vesicle->volume; 
    epsvol=V0*0.003e-2; //TODO! Follow Miha's derivation for exact formula;
    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++){
        vmsr=0.0;

 src/vertexmove.c

@@ -92,7 +92,7 @@
    memcpy((void *)&backupvtx[i+1],(void *)vtx->neigh[i],sizeof(ts_vertex));
    }

    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch==1){
    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch>0){
        for(i=0;i<vtx->tristar_no;i++) dvol-=vtx->tristar[i]->volume;
    };

@@ -113,11 +113,27 @@
        delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
    }

    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch == 1){
    if(vesicle->pswitch == 1 || vesicle->tape->constvolswitch >0){
        for(i=0;i<vtx->tristar_no;i++) dvol+=vtx->tristar[i]->volume;
        if(vesicle->pswitch==1) delta_energy-=vesicle->pressure*dvol;
    };


    if(vesicle->tape->constvolswitch==2){
        /*check whether the dvol is gt than epsvol */
            //fprintf(stderr,"DVOL=%1.16e\n",dvol);
        if(fabs(vesicle->volume+dvol-V0)>epsvol){
            //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;
        }

    } else
//    vesicle_volume(vesicle);
//    fprintf(stderr,"Volume before=%1.16e\n", vesicle->volume);
   if(vesicle->tape->constvolswitch == 1){
@@ -189,6 +205,9 @@
        
    }

    if(vesicle->tape->constvolswitch == 2){
    vesicle->volume+=dvol;
    } else
    if(vesicle->tape->constvolswitch == 1){
        constvolumeaccept(vesicle,constvol_vtx_moved,constvol_vtx_backup);
    }