| | |
| | | /* vim: set ts=4 sts=4 sw=4 noet : */ |
| | | #include<stdlib.h> |
| | | #include<math.h> |
| | | #include "general.h" |
| | |
| | | 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; |
| | |
| | | // 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; |
| | |
| | | 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; |
| | | } |
| | | } |
| | | |
| | |
| | | |
| | | // TODO: Maybe faster if checks only nucleus-neighboring cells |
| | | // Nucleus penetration check: |
| | | if(vesicle->R_nucleus>0.0){ |
| | | if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){ |
| | | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| | | return TS_FAIL; |
| | | } |
| | | } else if(vesicle->R_nucleusX>0.0){ |
| | | // fprintf(stderr,"DEBUG, (Rx, Ry,Rz)^2=(%f,%f,%f)\n",vesicle->R_nucleusX, vesicle->R_nucleusY, vesicle->R_nucleusZ); |
| | | if ((vtx->x*vtx->x)/vesicle->R_nucleusX + vtx->y*vtx->y/vesicle->R_nucleusY + (vtx->z*vtx->z)/vesicle->R_nucleusZ < 1.0){ |
| | | vtx=memcpy((void *)vtx,(void *)&backupvtx[0],sizeof(ts_vertex)); |
| | | return TS_FAIL; |
| | | } |
| | | |
| | | } |
| | | |
| | | //self avoidance check with distant vertices |
| | | cellidx=vertex_self_avoidance(vesicle, vtx); |
| | |
| | | } |
| | | |
| | | |
| | | //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; |
| | | |
| | |
| | | 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); |