Added direct force plugin and fixes in stretching energy

Samo Penic

2019-03-09 2c4278db6ead5c27e30a3000097ed898c968534e

 src/vertexmove.c

@@ -20,7 +20,7 @@
    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;
@@ -62,12 +62,6 @@
   }
/* 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]);
   oenergy=vtx->energy;
@@ -80,29 +74,16 @@
        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 +91,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 +134,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 */