From 17fe35ccc428e18dd226e07d5517c4816ef6be44 Mon Sep 17 00:00:00 2001
From: Samo Penic <samo.penic@gmail.com>
Date: Mon, 19 Apr 2021 12:23:07 +0000
Subject: [PATCH] Work done previously
---
src/tape | 6 +-
src/energy.c | 50 ++++++++++++++++++-------
src/initial_distribution.c | 4 ++
src/vertexmove.c | 1
4 files changed, 44 insertions(+), 17 deletions(-)
diff --git a/src/energy.c b/src/energy.c
index 1c50d84..5913863 100644
--- a/src/energy.c
+++ b/src/energy.c
@@ -9,6 +9,18 @@
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
+
+
+
+int cmpfunc(const void *x, const void *y)
+{
+ double diff= fabs(*(double*)x) - fabs(*(double*)y);
+ if(diff<0) return 1;
+ else return -1;
+}
+
+
+
/** @brief Wrapper that calculates energy of every vertex in vesicle
*
* Function calculated energy of every vertex in vesicle. It can be used in
@@ -116,9 +128,11 @@
ts_double We;
ts_double Av, We_Av;
+ ts_double eigenval[3];
+
gsl_matrix *gsl_Sv=gsl_matrix_alloc(3,3);
- gsl_vector_complex *Sv_eigen=gsl_vector_complex_alloc(3);
- gsl_eigen_nonsymm_workspace *workspace=gsl_eigen_nonsymm_alloc(3);
+ gsl_vector *Sv_eigen=gsl_vector_alloc(3);
+ gsl_eigen_symm_workspace *workspace=gsl_eigen_symm_alloc(3);
ts_double mprod[7], phi[7], he[7];
ts_double Sv[3][3]={{0,0,0},{0,0,0},{0,0,0}};
@@ -222,7 +236,7 @@
mprod[jj]=it->x*(k->y*edge_vector_z[jj]-edge_vector_y[jj]*k->z)-it->y*(k->x*edge_vector_z[jj]-k->z*edge_vector_x[jj])+it->z*(k->x*edge_vector_y[jj]-k->y*edge_vector_x[jj]);
- phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-8),mprod[jj])+M_PI;
+ phi[jj]=copysign(acos(lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm-1e-15),mprod[jj])+M_PI;
// printf("ACOS arg=%e\n", lm->xnorm*lp->xnorm+lm->ynorm*lp->ynorm+lm->znorm*lp->znorm);
//he was multiplied with 2 before...
he[jj]=sqrt( pow((edge_vector_x[jj]),2) + pow((edge_vector_y[jj]), 2) + pow((edge_vector_z[jj]), 2))*cos(phi[jj]/2.0);
@@ -265,21 +279,29 @@
// printf("Se= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Se11, Se21, Se31, Se21, Se22, Se32, Se31, Se32, Se33);
// printf("Pv= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Pv11, Pv21, Pv31, Pv21, Pv22, Pv32, Pv31, Pv32, Pv33);
- printf("Sv= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Sv[0][0], Sv[0][1], Sv[0][2], Sv[1][0], Sv[1][1], Sv[1][2], Sv[2][0], Sv[2][1], Sv[2][2]);
+// printf("Sv= %f, %f, %f\n %f, %f, %f\n %f, %f, %f\n", Sv[0][0], Sv[0][1], Sv[0][2], Sv[1][0], Sv[1][1], Sv[1][2], Sv[2][0], Sv[2][1], Sv[2][2]);
- gsl_eigen_nonsymm_params(0, 1, workspace);
- gsl_eigen_nonsymm(gsl_Sv, Sv_eigen, workspace);
- printf("Eigenvalues: %f, %f, %f\n",
- GSL_REAL(gsl_vector_complex_get(Sv_eigen, 0)),
- GSL_REAL(gsl_vector_complex_get(Sv_eigen, 1)),
- GSL_REAL(gsl_vector_complex_get(Sv_eigen, 2))
- );
- vtx->energy=0.0;
+ gsl_eigen_symm(gsl_Sv, Sv_eigen, workspace);
+
+// printf("Eigenvalues: %f, %f, %f\n", gsl_vector_get(Sv_eigen, 0),gsl_vector_get(Sv_eigen, 1), gsl_vector_get(Sv_eigen, 2) );
+// printf("Eigenvalues: %f, %f, %f\n", gsl_matrix_get(evec, 0,0),gsl_matrix_get(evec, 0,1), gsl_matrix_get(evec, 0,2) );
+
+
+ eigenval[0]= gsl_vector_get(Sv_eigen, 0);
+ eigenval[1]= gsl_vector_get(Sv_eigen, 1);
+ eigenval[2]= gsl_vector_get(Sv_eigen, 2);
+
+ qsort(eigenval, 3, sizeof(ts_double), cmpfunc);
+// printf("Eigenvalues: %f, %f, %f\n", eigenval[0], eigenval[1], eigenval[2] );
+
+
+ vtx->energy=(pow(eigenval[0]+eigenval[1],2))*Av;
gsl_matrix_free(gsl_Sv);
- gsl_vector_complex_free(Sv_eigen);
- gsl_eigen_nonsymm_free(workspace);
+ gsl_vector_free(Sv_eigen);
+// gsl_matrix_free(evec);
+ gsl_eigen_symm_free(workspace);
return TS_SUCCESS;
}
diff --git a/src/initial_distribution.c b/src/initial_distribution.c
index 87cfb68..bed8123 100644
--- a/src/initial_distribution.c
+++ b/src/initial_distribution.c
@@ -33,6 +33,10 @@
retval = mean_curvature_and_energy(vesicle);
ts_fprintf(stdout,"initial_distribution finished!\n");
if(retval);
+// ts_fprintf(stdout,"%e %e %e\n", vesicle->vlist->vtx[0]->x,vesicle->vlist->vtx[0]->y, vesicle->vlist->vtx[0]->z);
+// for(int i=0;i<vesicle->vlist->vtx[0]->neigh_no; i++){
+// ts_fprintf(stdout,"%e %e %e\n", vesicle->vlist->vtx[0]->neigh[i]->x,vesicle->vlist->vtx[0]->neigh[i]->y, vesicle->vlist->vtx[0]->neigh[i]->z);
+// }
return vesicle;
}
diff --git a/src/tape b/src/tape
index bdce9e7..4e78aff 100644
--- a/src/tape
+++ b/src/tape
@@ -14,7 +14,7 @@
# (pswitch=1: calc. p*dV energy contribution)
pswitch = 0
# pressure difference: p_inside - p_outside (in units kT/l_min^3):
-pressure=10.0
+pressure=10000.0
#Constant volume constraint (0 disable constant volume, 1 enable wiht additional vertex move, 2 enable with epsvol)
constvolswitch=0
@@ -30,7 +30,7 @@
####### Polymer (brush) definitions ###########
# npoly is a number of polymers attached to npoly distinct vertices on vesicle
-npoly=10
+npoly=0
# nmono is a number of monomers in each polymer
nmono=20
# Spring constant between monomers of the polymer
@@ -103,7 +103,7 @@
#plane confinement
-plane_confinement_switch=1
+plane_confinement_switch=0
#final plane distance (float in lmin)
plane_d=10
#plane to vesicle repulsion force while closing
diff --git a/src/vertexmove.c b/src/vertexmove.c
index e74fa74..c4ccb15 100644
--- a/src/vertexmove.c
+++ b/src/vertexmove.c
@@ -240,6 +240,7 @@
// fprintf(stderr, "DE=%f\n",delta_energy);
//MONTE CARLOOOOOOOO
// if(vtx->c!=0.0) printf("DE=%f\n",delta_energy);
+//delta_energy=1e15;
if(delta_energy>=0){
#ifdef TS_DOUBLE_DOUBLE
if(exp(-delta_energy)< drand48())
--
Gitblit v1.9.3