Trisurf Monte Carlo simulator
Samo Penic
2016-07-13 82a8abcddd62838bfcd06b6081042caf4d385a88
commit | author | age
7f6076 1 /* vim: set ts=4 sts=4 sw=4 noet : */
460c2a 2 #include<stdlib.h>
SP 3 #include<stdio.h>
4 #include<string.h>
5 #include<math.h>
6 #include "general.h"
7 #include "constvol.h"
8 #include "triangle.h"
9 #include "energy.h"
10 #include "vertex.h"
11 #include "cell.h"
12
fbcbdf 13 ts_bool constvolume(ts_vesicle *vesicle, ts_vertex *vtx_avoid, ts_double Vol, ts_double *retEnergy, ts_vertex **vtx_moved_retval, ts_vertex **vtx_backup){
SP 14     ts_vertex *vtx_moved;
460c2a 15     ts_uint vtxind,i,j;
41637a 16     ts_uint Ntries=100;
460c2a 17     ts_vertex *backupvtx;
fe5069 18     ts_double Rv, dh, dvol, volFirst, voldiff, oenergy,delta_energy;
460c2a 19     backupvtx=(ts_vertex *)calloc(sizeof(ts_vertex),10);
fbcbdf 20     ts_double l0 = (1.0 + sqrt(vesicle->dmax))/2.0; //make this a global constant if necessary
460c2a 21     for(i=0;i<Ntries;i++){
SP 22         vtxind=rand() % vesicle->vlist->n;
23         vtx_moved=vesicle->vlist->vtx[vtxind];
41637a 24         /* chosen vertex must not be a nearest neighbour. WASTODO: probably must.
SP 25          * DONE: solved below, when using different algorithm.
fe5069 26          * extend search in case of bondflip */
41637a 27 /*        if(vtx_moved==vtx_avoid) continue;
460c2a 28         for(j=0;j<vtx_moved->neigh_no;j++){
SP 29             if(vtx_moved->neigh[j]==vtx_avoid) continue;
30         }
41637a 31 */
SP 32 /* different check of nearest neighbour (and second nearest neighbour) check.
33  * Checking the distance between vertex and vertex to be moved to assure
34  * constant volume. Solves upper todo problem. */
35
36         if(vtx_distance_sq(vtx_moved,vtx_avoid)<16.0*vesicle->dmax){
37             continue;
38         }
39         
460c2a 40         memcpy((void *)&backupvtx[0],(void *)vtx_moved,sizeof(ts_vertex));
SP 41
fe5069 42         //move vertex in specified direction. first try, test move!
460c2a 43         Rv=sqrt(pow(vtx_moved->x,2)+pow(vtx_moved->y,2)+pow(vtx_moved->z,2));
fbcbdf 44         dh=2.0*Vol/(sqrt(3.0)*l0*l0);
SP 45         vtx_moved->x=vtx_moved->x*(1.0-dh/Rv);
46         vtx_moved->y=vtx_moved->y*(1.0-dh/Rv);
47         vtx_moved->z=vtx_moved->z*(1.0-dh/Rv);
460c2a 48
41637a 49
SP 50 //SKIPPING FIRST CHECK of CONSTRAINTS. This is not a final move.
460c2a 51         //check for constraints
41637a 52 /*          if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
460c2a 53             vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
41637a 54 //            continue;
SP 55                 break;
56         }  */
460c2a 57         // All checks OK!
SP 58
59         dvol=0.0;
60         for(j=0;j<vtx_moved->tristar_no;j++){
61             dvol-=vtx_moved->tristar[j]->volume;
fe5069 62         }
SP 63         volFirst=dvol;
64         for(j=0;j<vtx_moved->tristar_no;j++){
460c2a 65             triangle_normal_vector(vtx_moved->tristar[j]);
SP 66             dvol+=vtx_moved->tristar[j]->volume;
67         }
fe5069 68 //TODO: here there is a bug. Don't know where, but preliminary success can
SP 69 //happen sometimes. And when I do this checks, constant value is not achieved
70 //anymore
71 /*        voldiff=dvol-Vol;
460c2a 72         if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
SP 73             //calculate energy, return change in energy...
74              oenergy=vtx_moved->energy;
75             energy_vertex(vtx_moved);
76             delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
77             //the same is done for neighbouring vertices
fe5069 78             for(j=0;j<vtx_moved->neigh_no;j++){
SP 79                 oenergy=vtx_moved->neigh[j]->energy;
80                 energy_vertex(vtx_moved->neigh[j]);
81                 delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
460c2a 82             }
SP 83             *retEnergy=delta_energy;
fbcbdf 84             *vtx_backup=backupvtx;
SP 85             *vtx_moved_retval=vtx_moved;
b2fa8c 86             fprintf(stderr, "Preliminary success\n");
460c2a 87             return TS_SUCCESS;
fe5069 88         }        */
SP 89
90 //            fprintf(stderr, "Step 2 success\n");
460c2a 91         //do it again ;)
SP 92         dh=Vol*dh/dvol;
93         vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
94         vtx_moved->x=vtx_moved->x*(1-dh/Rv);
95         vtx_moved->y=vtx_moved->y*(1-dh/Rv);
96         vtx_moved->z=vtx_moved->z*(1-dh/Rv);
97         //check for constraints
98         if(constvolConstraintCheck(vesicle, vtx_moved)==TS_FAIL){
99             vtx_moved=memcpy((void *)vtx_moved,(void *)&backupvtx[0],sizeof(ts_vertex));
fe5069 100             //also, restore normals
SP 101             for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
41637a 102 //            continue;
SP 103                 break;
460c2a 104         }
SP 105
fe5069 106         dvol=volFirst;
b2fa8c 107         for(j=0;j<vtx_moved->tristar_no;j++){
SP 108             triangle_normal_vector(vtx_moved->tristar[j]);
109             dvol+=vtx_moved->tristar[j]->volume;
110         }
460c2a 111         voldiff=dvol-Vol;
SP 112         if(fabs(voldiff)/vesicle->volume < vesicle->tape->constvolprecision){
113             //calculate energy, return change in energy...
fe5069 114 //            fprintf(stderr, "Constvol success! %e\n",voldiff);
3de289 115             for(j=0;j<vtx_moved->neigh_no;j++){
SP 116             memcpy((void *)&backupvtx[j+1],(void *)vtx_moved->neigh[j],sizeof(ts_vertex));
117             }
118
460c2a 119             oenergy=vtx_moved->energy;
SP 120             energy_vertex(vtx_moved);
121             delta_energy=vtx_moved->xk*(vtx_moved->energy - oenergy);
122             //the same is done for neighbouring vertices
fe5069 123             for(j=0;j<vtx_moved->neigh_no;j++){
SP 124                 oenergy=vtx_moved->neigh[j]->energy;
125                 energy_vertex(vtx_moved->neigh[j]);
126                 delta_energy+=vtx_moved->neigh[j]->xk*(vtx_moved->neigh[j]->energy-oenergy);
460c2a 127             }
SP 128             *retEnergy=delta_energy;
fbcbdf 129             *vtx_backup=backupvtx;
SP 130             *vtx_moved_retval=vtx_moved;
460c2a 131             return TS_SUCCESS;
SP 132         }        
133     }
134     free(backupvtx);
135     return TS_FAIL;
136 }
137
138
139 ts_bool constvolConstraintCheck(ts_vesicle *vesicle, ts_vertex *vtx){ 
140         ts_uint i;
141         ts_double dist;
142         ts_uint cellidx;
143         //distance with neighbours check
144         for(i=0;i<vtx->neigh_no;i++){
145             dist=vtx_distance_sq(vtx,vtx->neigh[i]);
146             if(dist<1.0 || dist>vesicle->dmax) {
147             return TS_FAIL;
148             }
149         }
150         // Distance with grafted poly-vertex check:    
151         if(vtx->grafted_poly!=NULL){
152             dist=vtx_distance_sq(vtx,vtx->grafted_poly->vlist->vtx[0]);
153             if(dist<1.0 || dist>vesicle->dmax) {
154             return TS_FAIL;
155             }
156         }
157
158         // Nucleus penetration check:
159         if (vtx->x*vtx->x + vtx->y*vtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
160             return TS_FAIL;
161         }
162
163         //self avoidance check with distant vertices
164         cellidx=vertex_self_avoidance(vesicle, vtx);
165         //check occupation number
166         return cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
167 }
168
169
170
171 ts_bool constvolumerestore(ts_vertex *vtx_moved,ts_vertex *vtx_backup){
172     ts_uint j;
fbcbdf 173      memcpy((void *)vtx_moved,(void *)&vtx_backup[0],sizeof(ts_vertex));
fe5069 174     for(j=0;j<vtx_moved->tristar_no;j++) triangle_normal_vector(vtx_moved->tristar[j]);
3de289 175      for(j=0;j<vtx_moved->neigh_no;j++){
SP 176                // memcpy((void *)vtx_moved->neigh[j],(void *)&vtx_backup[j+1],sizeof(ts_vertex));
177                  energy_vertex(vtx_moved->neigh[j]);
178     }
fe5069 179
fbcbdf 180     free(vtx_backup);
460c2a 181     return TS_SUCCESS;
SP 182 }
183
fbcbdf 184 ts_bool constvolumeaccept(ts_vesicle *vesicle,ts_vertex *vtx_moved, ts_vertex *vtx_backup){
SP 185     ts_bool retval;
186     ts_uint cellidx=vertex_self_avoidance(vesicle, vtx_moved);
187     if(vtx_moved->cell!=vesicle->clist->cell[cellidx]){
188         retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx_moved);
189         if(retval==TS_SUCCESS) cell_remove_vertex(vtx_backup[0].cell,vtx_moved);
190         
191     }
192     free(vtx_backup);
460c2a 193
SP 194     return TS_SUCCESS;
195 }