From 3de2893d882d2c407375146376fbf9bb54f33e8b Mon Sep 17 00:00:00 2001 From: Samo Penic <samo.penic@fe.uni-lj.si> Date: Wed, 30 Apr 2014 13:44:12 +0000 Subject: [PATCH] Constant volume works. CAVEAT: in centermass I had to recalculate all normals of triangles again! This could influence behaviour of the system. --- src/constvol.c | 17 ++++++++--------- src/timestep.c | 10 ++++++++++ src/tape | 2 +- src/frame.c | 7 +++++++ src/vertexmove.c | 2 +- 5 files changed, 27 insertions(+), 11 deletions(-) diff --git a/src/constvol.c b/src/constvol.c index be69edf..72e902a 100644 --- a/src/constvol.c +++ b/src/constvol.c @@ -43,9 +43,6 @@ } // All checks OK! - for(j=0;j<vtx_moved->neigh_no;j++){ - memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex)); - } dvol=0.0; for(j=0;j<vtx_moved->tristar_no;j++){ dvol-=vtx_moved->tristar[j]->volume; @@ -86,9 +83,6 @@ vtx_moved->z=vtx_moved->z*(1-dh/Rv); //check for constraints if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){ - for(j=0;j<vtx_moved->neigh_no;j++){ - memcpy((void *)vtx_moved->neigh[j],(void *)&backupvtx[j+1],sizeof(ts_vertex)); - } vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex)); //also, restore normals for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]); @@ -104,6 +98,10 @@ if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){ //calculate energy, return change in energy... // fprintf(stderr, "Constvol success! %e\n",voldiff); + for(j=0;j<vtx_moved->neigh_no;j++){ + memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex)); + } + oenergy=vtx_moved->energy; energy_vertex(vtx_moved); delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy); @@ -159,10 +157,11 @@ ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){ ts_uint j; memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex)); - for(j=0;j<vtx_moved->neigh_no;j++){ - memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex)); - } for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]); + for(j=0;j<vtx_moved->neigh_no;j++){ + // memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex)); + energy_vertex(vtx_moved->neigh[j]); + } free(vtx_backup); return TS_SUCCESS; diff --git a/src/frame.c b/src/frame.c index 237aa87..ec94292 100644 --- a/src/frame.c +++ b/src/frame.c @@ -3,6 +3,8 @@ #include "cell.h" #include "frame.h" + +#include "triangle.h" ts_bool centermass(ts_vesicle *vesicle){ ts_uint i,j, n=vesicle->vlist->n; ts_vertex **vtx=vesicle->vlist->vtx; @@ -44,6 +46,11 @@ vesicle->cm[1]=0.0; vesicle->cm[2]=0.0; + for(i=0;i<vesicle->tlist->n;i++){ + triangle_normal_vector(vesicle->tlist->tria[i]); + } + + return TS_SUCCESS; } diff --git a/src/tape b/src/tape index 62074d5..ba80a8d 100644 --- a/src/tape +++ b/src/tape @@ -17,7 +17,7 @@ pressure=0.0 #Constant volume constraint (0 disable constant volume, 1 enable) -constvolswitch=1 +constvolswitch=0 constvolprecision=1e-14 diff --git a/src/timestep.c b/src/timestep.c index 7fd321c..a57687a 100644 --- a/src/timestep.c +++ b/src/timestep.c @@ -32,11 +32,16 @@ for(i=start_iteration;i<inititer+iterations;i++){ vmsr=0.0; bfsr=0.0; +/* vesicle_volume(vesicle); + fprintf(stderr,"Volume before TS=%1.16e\n", vesicle->volume); */ for(j=0;j<mcsweeps;j++){ single_timestep(vesicle, &vmsrt, &bfsrt); vmsr+=vmsrt; bfsr+=bfsrt; } +/* + vesicle_volume(vesicle); + fprintf(stderr,"Volume after TS=%1.16e\n", vesicle->volume); */ vmsr/=(ts_double)mcsweeps; bfsr/=(ts_double)mcsweeps; centermass(vesicle); @@ -87,6 +92,8 @@ } ts_bool single_timestep(ts_vesicle *vesicle,ts_double *vmsr, ts_double *bfsr){ +// vesicle_volume(vesicle); +// fprintf(stderr,"Volume before TS=%1.16e\n", vesicle->volume); ts_bool retval; ts_double rnvec[3]; ts_uint i,j, b; @@ -105,6 +112,7 @@ //find a bond and return a pointer to a bond... //call single_bondflip_timestep... retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec); + // b++; retval=TS_FAIL; if(retval==TS_SUCCESS) bfsrcnt++; } @@ -131,6 +139,8 @@ // printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,3*vesicle->blist->n,(double)cnt/(double)vesicle->blist->n/3.0); *vmsr=(ts_double)vmsrcnt/(ts_double)vesicle->vlist->n; *bfsr=(ts_double)bfsrcnt/(ts_double)vesicle->vlist->n/3.0; +// vesicle_volume(vesicle); +// fprintf(stderr,"Volume after TS=%1.16e\n", vesicle->volume); return TS_SUCCESS; } diff --git a/src/vertexmove.c b/src/vertexmove.c index cecb73e..af578e2 100644 --- a/src/vertexmove.c +++ b/src/vertexmove.c @@ -150,7 +150,7 @@ //MONTE CARLOOOOOOOO if(delta_energy>=0){ #ifdef TS_DOUBLE_DOUBLE - if(exp(-delta_energy)< drand48() ) + if(exp(-delta_energy)< drand48()) #endif #ifdef TS_DOUBLE_FLOAT if(expf(-delta_energy)< (ts_float)drand48()) -- Gitblit v1.9.3