From 42190fb471120cd6aafc9227081c57f7a0aca68e Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Thu, 04 Jul 2019 20:04:10 +0000
Subject: [PATCH] Fixed the max_occupancy double -> int

---
 src/poly.c |  159 +++++++++++++++++++++++++++++++++++++++++-----------
 1 files changed, 125 insertions(+), 34 deletions(-)

diff --git a/src/poly.c b/src/poly.c
index 3c93d0c..bcf01e5 100644
--- a/src/poly.c
+++ b/src/poly.c
@@ -6,7 +6,8 @@
 #include"bond.h"
 #include<math.h>
 #include"energy.h"
-
+#include"cell.h"
+#include"frame.h"
 ts_bool poly_assign_filament_xi(ts_vesicle *vesicle, ts_tape *tape){
 	ts_uint i;
 
@@ -47,8 +48,107 @@
 	poly->blist->bond[i]->bond_length=sqrt(vtx_distance_sq(poly->blist->bond[i]->vtx1,poly->blist->bond[i]->vtx2));
 	bond_energy(poly->blist->bond[i],poly);
 	}
-
+	vertex_list_assign_id(poly->vlist,TS_ID_FILAMENT);
 	return poly;
+}
+
+
+
+ts_bool poly_initial_distribution(ts_poly_list *poly_list, ts_int i, ts_vesicle *vesicle){
+	/* Make straight grafted poylmers normal to membrane (polymer brush). Dist. between poly vertices put to 1*/
+		ts_double xnorm,ynorm,znorm,normlength;
+		ts_int intpoly=vesicle->tape->internal_poly;
+		ts_int cellidx;
+		ts_double posX,posY,posZ,prevPosX,prevPosY,prevPosZ, phi,costheta,sintheta;
+		ts_bool retval;
+		ts_int j,k,l,m;
+			xnorm=0.0;
+			ynorm=0.0;
+			znorm=0.0;
+			for (j=0;j<poly_list->poly[i]->grafted_vtx->tristar_no;j++){
+				xnorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->xnorm;
+				ynorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->ynorm;
+				znorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->znorm;	
+			}
+			normlength=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
+			if(intpoly && i%2){
+				normlength=-normlength;
+			}
+			xnorm=xnorm/normlength;
+			ynorm=ynorm/normlength;
+			znorm=znorm/normlength;
+			//prepare starting position for building the polymeres
+			prevPosX=poly_list->poly[i]->grafted_vtx->x;
+			prevPosY=poly_list->poly[i]->grafted_vtx->y;
+			prevPosZ=poly_list->poly[i]->grafted_vtx->z;
+			for (j=0;j<poly_list->poly[i]->vlist->n;j++){
+				//if(j==0){
+					posX=prevPosX+xnorm*(vesicle->clist->dmin_interspecies);
+					posY=prevPosY+ynorm*(vesicle->clist->dmin_interspecies);
+					posZ=prevPosZ+znorm*(vesicle->clist->dmin_interspecies);
+				//}else{
+				//	posX=prevPosX+xnorm;
+				//	posY=prevPosY+ynorm;
+				//	posZ=prevPosZ+znorm;
+				//}
+				//trying to go towards normal
+				k=0;
+				while(1){
+					poly_list->poly[i]->vlist->vtx[j]->x = posX;
+					poly_list->poly[i]->vlist->vtx[j]->y = posY;
+					poly_list->poly[i]->vlist->vtx[j]->z = posZ;
+					cellidx=vertex_self_avoidance(vesicle, poly_list->poly[i]->vlist->vtx[j]);
+					retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,poly_list->poly[i]->vlist->vtx[j]);
+					if(retval==TS_SUCCESS){
+						retval=cell_add_vertex(vesicle->clist->cell[cellidx],poly_list->poly[i]->vlist->vtx[j]);
+						break;
+					}
+					else{
+					//	printf("%d %d Cannot put the vertex here. Finding another position\n",i,j);
+						//randomly change the direction.
+						m=0;
+						//we must move first vertex into the vesicle if the normal is in or out of the vesicle if the normal is out
+						do{
+							costheta=2.0*drand48()-1.0;
+							sintheta=sqrt(1-pow(costheta,2));
+							phi=drand48()*2.0*M_PI;
+							if(j==0){
+						//for special cases, when we are on the edge of bipyramid  the distance od dmin_interspecies is not enough
+								posX=prevPosX+vesicle->dmax*sintheta*cos(phi);
+								posY=prevPosY+vesicle->dmax*sintheta*sin(phi);
+								posZ=prevPosZ+vesicle->dmax*costheta;
+
+							} else {
+								posX=prevPosX+vesicle->clist->dmin_interspecies*sintheta*cos(phi);
+								posY=prevPosY+vesicle->clist->dmin_interspecies*sintheta*sin(phi);
+								posZ=prevPosZ+vesicle->clist->dmin_interspecies*costheta;
+							}
+							m++;
+							if(m>1000) {
+								k=9999; //break also ot of the outer loop
+								printf("was here\n");
+								break;
+							}
+						}
+						while((xnorm*(poly_list->poly[i]->grafted_vtx->x-posX)+ynorm*(poly_list->poly[i]->grafted_vtx->y-posY)+znorm*(poly_list->poly[i]->grafted_vtx->z-posZ))>0.0 && j==0);
+					}
+					k++;
+					if(k>1000){
+						//undo changes to the cell
+						for(l=0;l<j;l++){
+							cellidx=vertex_self_avoidance(vesicle, poly_list->poly[i]->vlist->vtx[l]);
+							cell_remove_vertex(vesicle->clist->cell[cellidx],poly_list->poly[i]->vlist->vtx[l]);
+						}
+						return TS_FAIL;
+					}
+				}
+				prevPosX=posX;
+				prevPosY=posY;
+				prevPosZ=posZ;
+			}
+	printf("did it\n");
+	return TS_SUCCESS;
+
 }
 
 
@@ -57,18 +157,30 @@
 	poly_list->poly	= (ts_poly **)calloc(n_poly,sizeof(ts_poly *));
 	ts_uint i=0,j=0; //idx;
 	ts_uint gvtxi;
-	ts_double xnorm,ynorm,znorm,normlength;
+	ts_bool retval;
 	ts_double dphi,dh;
-
+	cell_occupation(vesicle); //needed for evading the membrane
 	// Grafting polymers:
+	int tries=0;
 	if (vlist!=NULL){
 		if (n_poly > vlist->n){fatal("Number of polymers larger than numbero f vertices on a vesicle.",310);}
-	
 		while(i<n_poly){
 			gvtxi = rand() % vlist->n;
 			if (vlist->vtx[gvtxi]->grafted_poly == NULL){
-			poly_list->poly[i] = init_poly(n_mono, vlist->vtx[gvtxi]);
-			i++;
+				poly_list->poly[i] = init_poly(n_mono, vlist->vtx[gvtxi]);
+				retval=poly_initial_distribution(poly_list, i, vesicle);
+				if(retval==TS_FAIL){
+					ts_fprintf(stdout,"Found new potential grafting vertex %d for poly %d\n",gvtxi,i);
+					poly_free(poly_list->poly[i]);
+					tries++;
+				}
+				else {
+					tries=0;
+					i++;
+				}
+				if(tries>5000){
+					fatal("Cannot find space for inner polymeres",1001);
+				}
 			}
 		}
 	}
@@ -82,28 +194,6 @@
 	poly_list->n = n_poly;
 
 	if (vlist!=NULL){
-	/* Make straight grafted poylmers normal to membrane (polymer brush). Dist. between poly vertices put to 1*/
-		for (i=0;i<poly_list->n;i++){
-	
-			xnorm=0.0;
-			ynorm=0.0;
-			znorm=0.0;
-			for (j=0;j<poly_list->poly[i]->grafted_vtx->tristar_no;j++){
-				xnorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->xnorm;
-				ynorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->ynorm;
-				znorm-=poly_list->poly[i]->grafted_vtx->tristar[j]->znorm;	
-			}
-			normlength=sqrt(xnorm*xnorm+ynorm*ynorm+znorm*znorm);
-			xnorm=xnorm/normlength;
-			ynorm=ynorm/normlength;
-			znorm=znorm/normlength;
-
-			for (j=0;j<poly_list->poly[i]->vlist->n;j++){
-				poly_list->poly[i]->vlist->vtx[j]->x = poly_list->poly[i]->grafted_vtx->x + xnorm*(ts_double)(j+1);
-				poly_list->poly[i]->vlist->vtx[j]->y = poly_list->poly[i]->grafted_vtx->y + ynorm*(ts_double)(j+1);
-				poly_list->poly[i]->vlist->vtx[j]->z = poly_list->poly[i]->grafted_vtx->z + znorm*(ts_double)(j+1);
-			}
-		}
 	}
 	else
 	{
@@ -145,16 +235,16 @@
 		}
 	}
 */
-
+	
 	return poly_list;
 }
 
 
 ts_bool poly_free(ts_poly *poly){
 
-	if (poly->grafted_vtx!=NULL){
-		poly->grafted_vtx->grafted_poly=NULL;
-	}
+//	if (poly->grafted_vtx!=NULL){
+//		poly->grafted_vtx->grafted_poly=NULL;
+//	}
 	vtx_list_free(poly->vlist);
 	bond_list_free(poly->blist);
 	free(poly);
@@ -164,8 +254,9 @@
 
 ts_bool poly_list_free(ts_poly_list *poly_list){
 	ts_uint i;
-
+	//fprintf(stderr,"no. of polys=%d\n", poly_list->n);
 	for(i=0;i<poly_list->n;i++){
+	//	fprintf(stderr,"%d poly address in mem=%ld\n",i+1,(long)&(poly_list->poly[i]));
 		poly_free(poly_list->poly[i]);
 	}
 	free(poly_list->poly);

--
Gitblit v1.9.3