Merge branch 'master' of bitbucket.org:samop/trisurf-ng
| | |
| | | # generated automatically by aclocal 1.11.3 -*- Autoconf -*- |
| | | # generated automatically by aclocal 1.11.6 -*- Autoconf -*- |
| | | |
| | | # Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, |
| | | # 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, |
| | |
| | | |
| | | m4_ifndef([AC_AUTOCONF_VERSION], |
| | | [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl |
| | | m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],, |
| | | [m4_warning([this file was generated for autoconf 2.68. |
| | | m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.69],, |
| | | [m4_warning([this file was generated for autoconf 2.69. |
| | | You have another version of autoconf. It may work, but is not guaranteed to. |
| | | If you have problems, you may need to regenerate the build system entirely. |
| | | To do so, use the procedure documented by the package, typically `autoreconf'.])]) |
| | |
| | | [am__api_version='1.11' |
| | | dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to |
| | | dnl require some minimum version. Point them to the right macro. |
| | | m4_if([$1], [1.11.3], [], |
| | | m4_if([$1], [1.11.6], [], |
| | | [AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl |
| | | ]) |
| | | |
| | |
| | | # Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced. |
| | | # This function is AC_REQUIREd by AM_INIT_AUTOMAKE. |
| | | AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION], |
| | | [AM_AUTOMAKE_VERSION([1.11.3])dnl |
| | | [AM_AUTOMAKE_VERSION([1.11.6])dnl |
| | | m4_ifndef([AC_AUTOCONF_VERSION], |
| | | [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl |
| | | _AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))]) |
| | |
| | | //NOW insert vertices into data! |
| | | blist->bond[blist->n - 1]->vtx1=vtx1; |
| | | blist->bond[blist->n - 1]->vtx2=vtx2; |
| | | |
| | | blist->bond[blist->n - 1]->tainted=0; |
| | | //Should we calculate bond length NOW? |
| | | |
| | | return blist->bond[blist->n-1]; |
| | |
| | | ncy=(ts_uint)((vtx->y-vesicle->cm[1])*clist->dcell+clist->shift); |
| | | ncz=(ts_uint)((vtx->z-vesicle->cm[2])*clist->dcell+clist->shift); |
| | | |
| | | if(ncx == clist->ncmax[0]-1 || ncx == 2){ |
| | | if(ncx >= clist->ncmax[0]-1 || ncx <= 2){ |
| | | fatal("Vesicle is positioned outside the cell covered area. Coordinate x is the problem.",1500); |
| | | } |
| | | if(ncy == clist->ncmax[1]-1 || ncy == 2){ |
| | | if(ncy >= clist->ncmax[1]-1 || ncy <= 2){ |
| | | fatal("Vesicle is positioned outside the cell covered area. Coordinate y is the problem.",1500); |
| | | } |
| | | if(ncz == clist->ncmax[2]-1 || ncz == 2){ |
| | | if(ncz >= clist->ncmax[2]-1 || ncz <= 2){ |
| | | fatal("Vesicle is positioned outside the cell covered area. Coordinate z is the problem.",1500); |
| | | } |
| | | cellidx=ncz+(ncy-1)*clist->ncmax[2] + (ncx-1)*clist->ncmax[2]* |
| | |
| | | |
| | | |
| | | //TODO: looks ok, but debug anyway in the future |
| | | ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){ |
| | | inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx){ |
| | | ts_uint i; |
| | | for(i=0;i<cell->nvertex;i++){ |
| | | if(cell->vertex[i]==vtx){ |
| | | //vertex is already in the cell! |
| | | //fprintf(stderr,"VTX in the cell!\n"); |
| | | return TS_FAIL; |
| | | } |
| | | } |
| | | //fprintf(stderr,"VTX added to the cell!\n"); |
| | | cell->nvertex++; |
| | | cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *)); |
| | | if(cell->vertex == NULL){ |
| | | fatal("Reallocation of memory failed during insertion of vertex in cell_add_vertex",3); |
| | | } |
| | | cell->vertex[cell->nvertex-1]=vtx; |
| | | vtx->cell=cell; |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | | inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx){ |
| | | ts_uint i,j=0; |
| | | for(i=0;i<cell->nvertex;i++){ |
| | | if(cell->vertex[i]!=vtx){ |
| | | cell->vertex[j]=cell->vertex[i]; |
| | | j++; |
| | | } |
| | | } |
| | | if(j==i){ |
| | | fatal("Vertex was not in the cell!",3); |
| | | } |
| | | //fprintf(stderr, "Vertex deleted from the cell!\n"); |
| | | |
| | | /* resize memory. potentionally time consuming */ |
| | | cell->nvertex--; |
| | | cell->vertex=(ts_vertex **)realloc(cell->vertex,cell->nvertex*sizeof(ts_vertex *)); |
| | | if(vtx->neigh == NULL && vtx->neigh_no!=0) |
| | | if(cell->vertex == NULL){ |
| | | fatal("Reallocation of memory failed during removal of vertex in cell_remove_vertex",3); |
| | | } |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | | ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist){ |
| | | ts_uint i; |
| | |
| | | } |
| | | |
| | | // TODO: compiles ok, but it is completely untested and undebugged. It was debugged before rewrite, but this was long time ago. |
| | | ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx){ |
| | | ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx){ |
| | | ts_uint ncx,ncy,ncz,remainder,cell_occupation; |
| | | ts_uint i,j,k,l,neigh_cidx; |
| | | ts_double dist; |
| | |
| | | // cell! |
| | | if(cell_occupation>1){ |
| | | for(l=0;l<cell_occupation;l++){ |
| | | if(clist->cell[neigh_cidx]->vertex[l]!=vtx){ |
| | | // fprintf(stderr,"calling dist on vertex %i\n",l); |
| | | dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],tvtx); |
| | | // fprintf(stderr,"dist was %f\n",dist); |
| | | if(dist<1) return TS_FAIL; |
| | | } |
| | | } |
| | | } |
| | | } |
| | | } |
| | | } |
| | | |
| | | //carefull with this checks! |
| | | if(clist->cell[neigh_cidx]->vertex[l]->idx!=vtx->idx){ |
| | | // fprintf(stderr,"calling dist on vertex %i\n",l); |
| | | dist=vtx_distance_sq(clist->cell[neigh_cidx]->vertex[l],vtx); |
| | | // fprintf(stderr,"dist was %f\n",dist); |
| | | if(dist<=1.0) return TS_FAIL; |
| | | } |
| | | } |
| | | } |
| | | } |
| | | } |
| | | } |
| | | return TS_SUCCESS; |
| | | } |
| | |
| | | ts_bool cell_free(ts_cell* cell); |
| | | ts_bool cell_list_free(ts_cell_list *clist); |
| | | inline ts_uint vertex_self_avoidance(ts_vesicle *vesicle, ts_vertex *vtx); |
| | | ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx); |
| | | inline ts_bool cell_add_vertex(ts_cell *cell, ts_vertex *vtx); |
| | | inline ts_bool cell_remove_vertex(ts_cell *cell, ts_vertex *vtx); |
| | | ts_bool cell_list_cell_occupation_clear(ts_cell_list *clist); |
| | | |
| | | ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx, ts_vertex *tvtx); |
| | | inline ts_bool cell_occupation_number_and_internal_proximity(ts_cell_list *clist, ts_uint cellidx, ts_vertex *vtx); |
| | | #endif |
| | |
| | | ts_uint jjp,jjm; |
| | | ts_vertex *j,*jp, *jm; |
| | | ts_triangle *jt; |
| | | ts_double s=0,xh=0,yh=0,zh=0,txn=0,tyn=0,tzn=0; |
| | | ts_double s=0.0,xh=0.0,yh=0.0,zh=0.0,txn=0.0,tyn=0.0,tzn=0.0; |
| | | ts_double x1,x2,x3,ctp,ctm,tot,xlen; |
| | | ts_double h,ht; |
| | | for(jj=1; jj<=vtx->neigh_no;jj++){ |
| | |
| | | #endif |
| | | tot=ctp+ctm; |
| | | tot=0.5*tot; |
| | | |
| | | xlen=vtx_distance_sq(j,vtx); |
| | | /* |
| | | #ifdef TS_DOUBLE_DOUBLE |
| | | vtx->bond[jj-1]->bond_length=sqrt(xlen); |
| | | #endif |
| | |
| | | #endif |
| | | |
| | | vtx->bond[jj-1]->bond_length_dual=tot*vtx->bond[jj-1]->bond_length; |
| | | |
| | | */ |
| | | s+=tot*xlen; |
| | | xh+=tot*(j->x - vtx->x); |
| | | yh+=tot*(j->y - vtx->y); |
| | |
| | | vtx[i]->z-=vesicle->cm[2]; |
| | | } |
| | | |
| | | vesicle->cm[0]=0; |
| | | vesicle->cm[1]=0; |
| | | vesicle->cm[2]=0; |
| | | |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | |
| | | cell_list_cell_occupation_clear(vesicle->clist); |
| | | for(i=0;i<n;i++){ |
| | | cellidx=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]); |
| | | vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx]; |
| | | // already done in cell_add_vertex |
| | | // vesicle->vlist->vtx[i]->cell=vesicle->clist->cell[cellidx]; |
| | | |
| | | cell_add_vertex(vesicle->clist->cell[cellidx],vesicle->vlist->vtx[i]); |
| | | |
| | |
| | | ts_vertex *vtx2; |
| | | ts_double bond_length; |
| | | ts_double bond_length_dual; |
| | | ts_bool tainted; |
| | | }; |
| | | typedef struct ts_bond ts_bond; |
| | | |
| | |
| | | #include<stdlib.h> |
| | | #include <sys/types.h> |
| | | #include <dirent.h> |
| | | |
| | | #include "initial_distribution.h" |
| | | |
| | | ts_bool print_vertex_list(ts_vertex_list *vlist){ |
| | | ts_uint i; |
| | |
| | | |
| | | |
| | | |
| | | ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations){ |
| | | ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations){ |
| | | long int nshell=17,ncxmax=60, ncymax=60, nczmax=60; // THIS IS DUE TO CONFUSE BUG! |
| | | char buf[255]; |
| | | long int brezveze=1; |
| | | char *buf=malloc(255*sizeof(char)); |
| | | long int brezveze0=1; |
| | | long int brezveze1=1; |
| | | long int brezveze2=1; |
| | | ts_double xk0=25.0, dmax=1.67,stepsize=0.15; |
| | | *iterations=1000; |
| | | long int iter=1000, init=1000, mcsw=1000; |
| | | cfg_opt_t opts[] = { |
| | | CFG_SIMPLE_INT("nshell", &nshell), |
| | | CFG_SIMPLE_FLOAT("dmax", &dmax), |
| | |
| | | CFG_SIMPLE_INT("nxmax", &ncxmax), |
| | | CFG_SIMPLE_INT("nymax", &ncymax), |
| | | CFG_SIMPLE_INT("nzmax", &nczmax), |
| | | CFG_SIMPLE_INT("iterations",iterations), |
| | | CFG_SIMPLE_INT("iterations",&iter), |
| | | CFG_SIMPLE_INT("mcsweeps",&mcsw), |
| | | CFG_SIMPLE_INT("inititer", &init), |
| | | CFG_SIMPLE_BOOL("quiet",&quiet), |
| | | CFG_SIMPLE_STR("multiprocessing",buf), |
| | | CFG_SIMPLE_INT("smp_cores",&brezveze), |
| | | CFG_SIMPLE_INT("cluster_nodes",&brezveze), |
| | | CFG_SIMPLE_INT("distributed_processes",&brezveze), |
| | | CFG_SIMPLE_INT("smp_cores",&brezveze0), |
| | | CFG_SIMPLE_INT("cluster_nodes",&brezveze1), |
| | | CFG_SIMPLE_INT("distributed_processes",&brezveze2), |
| | | CFG_END() |
| | | }; |
| | | cfg_t *cfg; |
| | |
| | | else if(retval==CFG_PARSE_ERROR){ |
| | | fatal("Invalid tape!",100); |
| | | } |
| | | ts_vesicle *vesicle; |
| | | *iterations=iter; |
| | | *inititer=init; |
| | | *mcsweeps=mcsw; |
| | | vesicle=initial_distribution_dipyramid(nshell,ncxmax,ncymax,nczmax,stepsize); |
| | | vesicle->nshell=nshell; |
| | | vesicle->dmax=dmax*dmax; |
| | | vesicle->bending_rigidity=xk0; |
| | |
| | | vesicle->clist->ncmax[1]=ncymax; |
| | | vesicle->clist->ncmax[2]=nczmax; |
| | | vesicle->clist->max_occupancy=8; |
| | | |
| | | cfg_free(cfg); |
| | | free(buf); |
| | | // fprintf(stderr,"NSHELL=%u\n",vesicle->nshell); |
| | | return TS_SUCCESS; |
| | | return vesicle; |
| | | |
| | | } |
| | | |
| | |
| | | ts_bool write_vertex_vtk_file(ts_vesicle *vesicle,ts_char *filename, ts_char *text); |
| | | ts_bool write_vertex_xml_file(ts_vesicle *vesicle, ts_uint timestepno); |
| | | ts_bool write_master_xml_file(ts_char *filename); |
| | | ts_bool parsetape(ts_vesicle *vesicle,ts_uint *iterations); |
| | | |
| | | ts_vesicle *parsetape(ts_uint *mcsweeps, ts_uint *inititer, ts_uint *iterations); |
| | | |
| | | #endif |
| | |
| | | */ |
| | | |
| | | int main(int argv, char *argc[]){ |
| | | ts_uint inititer,mcsweeps, iterations; |
| | | ts_vesicle *vesicle; |
| | | /* THIS SHOULD GO INTO UNIT TEST |
| | | ts_bool retval; |
| | | ts_uint i; |
| | | ts_vertex_list *vlist=init_vertex_list(5); |
| | | ts_vertex_list *vlist1; |
| | | ts_bond_list *blist=init_bond_list(); |
| | | ts_triangle_list *tlist=init_triangle_list(); |
| | | ts_cell_list *clist=init_cell_list(3,3,3,0.3); |
| | | ts_vesicle *vesicle; |
| | | |
| | | retval=vtx_add_cneighbour(blist,vlist->vtx[1],vlist->vtx[0]); |
| | | if(retval==TS_FAIL) printf("1. already a member or vertex is null!\n"); |
| | |
| | | |
| | | vtx_list_free(vlist1); |
| | | printf("Tests complete.\n"); |
| | | vesicle=initial_distribution_dipyramid(17,60,60,60,0.15); |
| | | //parsetape(vesicle,&i); |
| | | */ |
| | | vesicle=parsetape(&mcsweeps, &inititer, &iterations); |
| | | run_simulation(vesicle, mcsweeps, inititer, iterations); |
| | | |
| | | //these four must come from parsetype! |
| | | vesicle->dmax=1.67*1.67; |
| | | vesicle->stepsize=0.15; |
| | | vesicle->clist->max_occupancy=8; |
| | | vesicle->bending_rigidity=25.0; |
| | | fprintf(stderr,"xk=%f",vesicle->bending_rigidity); |
| | | |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | for(i=0;i<100;i++){ |
| | | single_timestep(vesicle); |
| | | if(i%100==0){ |
| | | write_vertex_xml_file(vesicle,i/100); |
| | | } |
| | | } |
| | | write_master_xml_file("test.pvd"); |
| | | write_dout_fcompat_file(vesicle,"dout"); |
| | | vesicle_free(vesicle); |
| | |
| | | */ |
| | | ts_bool saveAvgUlm2(ts_vesicle *vesicle); |
| | | int main(int argv, char *argc[]){ |
| | | ts_uint i,j; |
| | | ts_uint i,j,k; |
| | | ts_vesicle *vesicle; |
| | | ts_double r0; |
| | | vesicle=initial_distribution_dipyramid(17,60,60,60,0.15); |
| | | //parsetape(vesicle,&i); |
| | | |
| | | //similar to nmax in fortran code |
| | | ts_uint nmax; |
| | | |
| | | //these four must come from parsetype! |
| | | vesicle->dmax=1.67*1.67; |
| | |
| | | //fprintf(stderr,"xk=%f",vesicle->bending_rigidity); |
| | | |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | |
| | | //test if the structure is internally organized into cells correctly |
| | | ts_uint cind; |
| | | for(i=0;i<vesicle->vlist->n;i++){ |
| | | cind=vertex_self_avoidance(vesicle, vesicle->vlist->vtx[i]); |
| | | |
| | | if(vesicle->clist->cell[cind]==vesicle->vlist->vtx[i]->cell){ |
| | | //fprintf(stdout,"(T) Idx match!\n"); |
| | | } else { |
| | | fprintf(stderr,"(T) ***** Idx don't match!\n"); |
| | | |
| | | } |
| | | } |
| | | //end test |
| | | vesicle->sphHarmonics=sph_init(vesicle->vlist, 21); |
| | | |
| | | vesicle_volume(vesicle); |
| | |
| | | calculateYlmi(vesicle); |
| | | calculateUlm(vesicle); |
| | | |
| | | |
| | | |
| | | for(i=0;i<10;i++){ |
| | | //preloop: |
| | | for(i=0;i<1000;i++){ |
| | | cell_occupation(vesicle); |
| | | for(j=0;j<1000;j++){ |
| | | single_timestep(vesicle); |
| | | } |
| | | centermass(vesicle); |
| | | fprintf(stderr, "Preloop %d completed.\n",i+1); |
| | | } |
| | | |
| | | nmax=1000; |
| | | for(i=0;i<nmax;i++){ |
| | | for(j=0;j<200;j++){ |
| | | cell_occupation(vesicle); |
| | | for(k=0;k<5;k++){ |
| | | single_timestep(vesicle); |
| | | } |
| | | centermass(vesicle); |
| | | } |
| | | vesicle_volume(vesicle); |
| | | r0=getR0(vesicle); |
| | | |
| | |
| | | saveAvgUlm2(vesicle); |
| | | |
| | | write_vertex_xml_file(vesicle,i); |
| | | fprintf(stderr, "Loop %d completed.\n",i+1); |
| | | fprintf(stderr, "Loop %d out of %d completed.\n",i+1,nmax); |
| | | |
| | | } |
| | | |
| | | write_master_xml_file("test.pvd"); |
| | | write_dout_fcompat_file(vesicle,"dout"); |
| | | vesicle_free(vesicle); |
| | |
| | | |
| | | |
| | | ####### Program Control ############ |
| | | #how many iteration are there in a run? |
| | | iterations=20000 |
| | | #how many MC sweeps between subsequent records of states to disk |
| | | mcsweeps=100 |
| | | #how many initial mcsweeps*inititer MC sweeps before recording to disk? |
| | | inititer=100 |
| | | #how many records do you want on the disk iteration are there in a run? |
| | | iterations=1000 |
| | | |
| | | |
| | | #shut up if we are using cluster!!! |
| | | quiet=false |
| | |
| | | #include "timestep.h" |
| | | #include "vertexmove.h" |
| | | #include "bondflip.h" |
| | | #include "frame.h" |
| | | #include "io.h" |
| | | ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){ |
| | | ts_uint i, j; |
| | | |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps); |
| | | for(i=0;i<inititer+iterations;i++){ |
| | | for(j=0;j<mcsweeps;j++){ |
| | | single_timestep(vesicle); |
| | | } |
| | | centermass(vesicle); |
| | | cell_occupation(vesicle); |
| | | if(i>inititer){ |
| | | write_vertex_xml_file(vesicle,i-inititer); |
| | | } |
| | | } |
| | | return TS_SUCCESS; |
| | | } |
| | | |
| | | ts_bool single_timestep(ts_vesicle *vesicle){ |
| | | ts_bool retval; |
| | | ts_double rnvec[3]; |
| | | ts_uint i; |
| | | ts_uint i, b; |
| | | for(i=0;i<vesicle->vlist->n;i++){ |
| | | rnvec[0]=drand48(); |
| | | rnvec[1]=drand48(); |
| | |
| | | } |
| | | |
| | | // ts_int cnt=0; |
| | | for(i=0;i<vesicle->blist->n;i++){ |
| | | rnvec[0]=drand48(); |
| | | for(i=0;i<vesicle->vlist->n;i++){ |
| | | //why is rnvec needed in bondflip? |
| | | /* rnvec[0]=drand48(); |
| | | rnvec[1]=drand48(); |
| | | rnvec[2]=drand48(); |
| | | */ |
| | | b=rand() % vesicle->blist->n; |
| | | //find a bond and return a pointer to a bond... |
| | | //call single_bondflip_timestep... |
| | | retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[i],rnvec); |
| | | retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec); |
| | | // if(retval==TS_SUCCESS) cnt++; |
| | | } |
| | | // printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n); |
| | |
| | | #ifndef _TIMESTEP_H |
| | | #define _TIMESTEP_H |
| | | ts_bool single_timestep(ts_vesicle *vesicle); |
| | | ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations); |
| | | #endif |
| | |
| | | //#include "io.h" |
| | | #include<stdio.h> |
| | | #include "vertexmove.h" |
| | | #include <string.h> |
| | | |
| | | ts_bool single_verticle_timestep(ts_vesicle *vesicle,ts_vertex *vtx,ts_double |
| | | *rn){ |
| | |
| | | ts_double dist; |
| | | ts_bool retval; |
| | | ts_uint cellidx; |
| | | ts_double xold,yold,zold; |
| | | ts_double delta_energy,oenergy; |
| | | ts_vertex *ovtx; |
| | | ts_vertex *tvtx=(ts_vertex *)calloc(1,sizeof(ts_vertex)); |
| | | ts_double costheta,sintheta,phi,r; |
| | | //This will hold all the information of vtx and its neighbours |
| | | ts_vertex backupvtx[20]; |
| | | memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex)); |
| | | |
| | | //randomly we move the temporary vertex |
| | | tvtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0); |
| | | tvtx->y=vtx->y+vesicle->stepsize*(2.0*rn[1]-1.0); |
| | | tvtx->z=vtx->z+vesicle->stepsize*(2.0*rn[2]-1.0); |
| | | //check we if some length to neighbours are too much |
| | | //Some stupid tests for debugging cell occupation! |
| | | /* cellidx=vertex_self_avoidance(vesicle, vtx); |
| | | if(vesicle->clist->cell[cellidx]==vtx->cell){ |
| | | fprintf(stderr,"Idx match!\n"); |
| | | } else { |
| | | fprintf(stderr,"***** Idx don't match!\n"); |
| | | fatal("ENding.",1); |
| | | } |
| | | */ |
| | | |
| | | //temporarly moving the vertex |
| | | // vtx->x=vtx->x+vesicle->stepsize*(2.0*rn[0]-1.0); |
| | | // 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: |
| | | r=vesicle->stepsize*rn[0]; |
| | | phi=rn[1]*2*M_PI; |
| | | costheta=2*rn[2]-1; |
| | | sintheta=sqrt(1-pow(costheta,2)); |
| | | vtx->x=vtx->x+r*sintheta*cos(phi); |
| | | vtx->y=vtx->y+r*sintheta*sin(phi); |
| | | vtx->z=vtx->z+r*costheta; |
| | | |
| | | |
| | | //distance with neighbours check |
| | | for(i=0;i<vtx->neigh_no;i++){ |
| | | dist=vtx_distance_sq(tvtx,vtx->neigh[i]); |
| | | dist=vtx_distance_sq(vtx,vtx->neigh[i]); |
| | | if(dist<1.0 || dist>vesicle->dmax) { |
| | | vtx_free(tvtx); |
| | | // fprintf(stderr,"Fail 1, dist=%f, vesicle->dmax=%f\n", dist, vesicle->dmax); |
| | | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| | | return TS_FAIL; |
| | | } |
| | | } |
| | | //self avoidance check with distant vertices |
| | | cellidx=vertex_self_avoidance(vesicle, tvtx); |
| | | cellidx=vertex_self_avoidance(vesicle, vtx); |
| | | //check occupation number |
| | | retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx,tvtx); |
| | | retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx); |
| | | |
| | | if(retval==TS_FAIL){ |
| | | vtx_free(tvtx); |
| | | // fprintf(stderr,"Fail 2\n"); |
| | | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| | | return TS_FAIL; |
| | | } |
| | | |
| | | //if all the tests are successful, then we update the vertex position |
| | | xold=vtx->x; |
| | | yold=vtx->y; |
| | | zold=vtx->z; |
| | | ovtx=malloc(sizeof(ts_vertex)); |
| | | vtx_copy(ovtx,vtx); |
| | | vtx->x=tvtx->x; |
| | | vtx->y=tvtx->y; |
| | | vtx->z=tvtx->z; |
| | | |
| | | //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)); |
| | | } |
| | | |
| | | |
| | | |
| | | delta_energy=0; |
| | | //update the normals of triangles that share bead i. |
| | | for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); |
| | | //energy and curvature |
| | | oenergy=vtx->energy; |
| | | energy_vertex(vtx); |
| | | delta_energy=vtx->xk*(vtx->energy - ovtx->energy); |
| | | delta_energy=vtx->xk*(vtx->energy - oenergy); |
| | | //the same is done for neighbouring vertices |
| | | for(i=0;i<vtx->neigh_no;i++){ |
| | | oenergy=vtx->neigh[i]->energy; |
| | |
| | | #endif |
| | | { |
| | | //not accepted, reverting changes |
| | | vtx->x=xold; |
| | | vtx->y=yold; |
| | | vtx->z=zold; |
| | | 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)); |
| | | } |
| | | |
| | | //update the normals of triangles that share bead i. |
| | | for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]); |
| | | //energy and curvature |
| | | energy_vertex(vtx); |
| | | //the same is done for neighbouring vertices |
| | | for(i=0;i<vtx->neigh_no;i++) energy_vertex(vtx->neigh[i]); |
| | | free(ovtx->bond_length); |
| | | free(ovtx->bond_length_dual); |
| | | free(ovtx); |
| | | vtx_free(tvtx); |
| | | |
| | | return TS_FAIL; |
| | | } |
| | | } |
| | | //END MONTE CARLOOOOOOO |
| | | |
| | | //TODO: change cell occupation if necessary! |
| | | // fprintf(stderr,"Success!!\n"); |
| | | free(ovtx->bond_length); |
| | | free(ovtx->bond_length_dual); |
| | | free(ovtx); |
| | | vtx_free(tvtx); |
| | | // oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]); |
| | | if(vtx->cell!=vesicle->clist->cell[cellidx]){ |
| | | retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx); |
| | | // if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx); |
| | | if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx); |
| | | |
| | | } |
| | | // if(oldcellidx); |
| | | //END MONTE CARLOOOOOOO |
| | | return TS_SUCCESS; |
| | | } |
| | | |