Trisurf Monte Carlo simulator
Samo Penic
2019-03-09 9a1f16ae3affc4db83f2eb2623418ac5cff2af03
src/vertexmove.c
@@ -20,10 +20,9 @@
    ts_uint i;
    ts_bool retval; 
    ts_uint cellidx; 
    ts_double delta_energy, oenergy, darea=0.0, dstretchenergy=0.0;
    ts_double delta_energy, oenergy;
    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;
    ts_vertex backupvtx[20];
    memcpy((void *)&backupvtx[0],(void *)vtx,sizeof(ts_vertex));
    //random move in a sphere with radius stepsize:
@@ -54,7 +53,6 @@
   }
/* Entry point for plugin vm_energy_before_prepare() */
   vesicle->plist->pointer=vesicle->plist->chain->vm_energy_before_prepare;
   while(vesicle->plist->pointer!=NULL){
      vesicle->plist->pointer->plugin->function->vm_energy_before_prepare(vesicle, vtx);
@@ -62,65 +60,27 @@
   }
/* End of vm_energy_before_prepare() */
   //stretching energy 1 of 3
   if(vesicle->tape->stretchswitch==1){
      for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
   }
    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]);
   //update the normals of triangles that share bead i.
   for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
   oenergy=vtx->energy;
    energy_vertex(vtx);
    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;
        energy_vertex(vtx->neigh[i]);
        delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
    }
   energy_vertex(vtx);
   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;
      energy_vertex(vtx->neigh[i]);
      delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
   }
/* Entry point for plugin vm_energy_after_execute() */
   vesicle->plist->pointer=vesicle->plist->chain->vm_energy_after_execute;
   while(vesicle->plist->pointer!=NULL){
      delta_energy+=vesicle->plist->pointer->plugin->function->vm_energy_after_execute(vesicle, vtx);
      delta_energy+=vesicle->plist->pointer->plugin->function->vm_energy_after_execute(vesicle, vtx, backupvtx);
      vesicle->plist->pointer=vesicle->plist->pointer->next;
   }
/* End of vm_energy_after_execute() */
    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;
      }
    }
/* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */
   delta_energy+=direct_force_energy(vesicle,vtx,backupvtx);
   //stretching energy 2 of 3
   if(vesicle->tape->stretchswitch==1){
      for(i=0;i<vtx->tristar_no;i++){
         stretchenergy(vesicle, vtx->tristar[i]);
         dstretchenergy+=vtx->tristar[i]->energy;
         }
   }
   delta_energy+=dstretchenergy;
/* No poly-bond energy for now!
   if(vtx->grafted_poly!=NULL){
      delta_energy+=
@@ -128,21 +88,6 @@
         pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
   }
*/
// plane confinement energy due to compressing force
   if(vesicle->tape->plane_confinement_switch){
      if(vesicle->confinement_plane.force_switch){
         //substract old energy
         if(abs(vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_max)>1e-10) {
            delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-backupvtx[0].z,2);
            delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-vtx->z,2);
         }
         if(abs(-vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_min)>1e-10) {
            delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-backupvtx[0].z,2);
            delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-vtx->z,2);
         }
      }
   }
/* Entry point for plugin vm_before_montecarlo_constraint() function */
@@ -187,14 +132,6 @@
   //update the normals of triangles that share bead i.
   for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
   //stretching energy 3 of 3
   if(vesicle->tape->stretchswitch==1){
      for(i=0;i<vtx->tristar_no;i++){
         stretchenergy(vesicle,vtx->tristar[i]);
         }
   }
/* Entry point for plugin vm_before_montecarlo_constraint() function */
   vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_rejected;
@@ -214,12 +151,6 @@
      retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
      if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);   
   }
    if(vesicle->tape->constareaswitch==2){
        vesicle->area+=darea;
    }
/* Entry point for plugin vm_before_montecarlo_constraint() function */
   vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_accepted;