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,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,47 +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() */
/* 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+=
@@ -110,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 */
@@ -168,14 +131,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 */