Trisurf Monte Carlo simulator
Samo Penic
2014-04-21 459ff94855b2659de14ac2ac902b9658d28e19d5
Added complex spherical harmonics.
2 files added
7 files modified
169 ■■■■■ changed files
src/Makefile.am 10 ●●●● patch | view | raw | blame | history
src/general.h 2 ●●●●● patch | view | raw | blame | history
src/initial_distribution.c 3 ●●●● patch | view | raw | blame | history
src/main.c 3 ●●●● patch | view | raw | blame | history
src/shcomplex.c 127 ●●●●● patch | view | raw | blame | history
src/shcomplex.h 8 ●●●●● patch | view | raw | blame | history
src/tape 6 ●●●● patch | view | raw | blame | history
src/timestep.c 7 ●●●●● patch | view | raw | blame | history
src/vesicle.c 3 ●●●● patch | view | raw | blame | history
src/Makefile.am
@@ -1,17 +1,17 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c bondflip.c main.c poly.c stats.c sh.c shcomplex.c
#trisurf_LDFLAGS = -lm -lconfuse
shdiscoverdir=../
shdiscover_PROGRAMS= shdiscover
shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c poly.c stats.c
shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c poly.c stats.c shcomplex.c
AM_CFLAGS = -Wall -Werror
co_testdir=../
co_test_PROGRAMS=co_test
co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c poly.c stats.c
co_test_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c co_test.c frame.c bondflip.c poly.c stats.c shcomplex.c
spherical_trisurfdir=../
spherical_trisurf_PROGRAMS = spherical_trisurf
spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c poly.c stats.c
spherical_trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf.c sh.c bondflip.c poly.c stats.c shcomplex.c
spherical_trisurf_ffdir=../
spherical_trisurf_ff_PROGRAMS = spherical_trisurf_ff
spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c poly.c stats.c
spherical_trisurf_ff_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c spherical_trisurf_ff.c sh.c bondflip.c poly.c stats.c shcomplex.c
src/general.h
@@ -3,6 +3,7 @@
#include<stdarg.h>
#include<stdio.h>
#include<gsl/gsl_complex.h>
/* @brief This is a header file, defining general constants and structures.
  * @file header.h
  * @author Samo Penic
@@ -222,6 +223,7 @@
typedef struct {
    ts_uint l;
    ts_double **ulm;
    gsl_complex **ulmComplex;
    ts_double **sumUlm2;
    ts_uint N;
    ts_double **co;
src/initial_distribution.c
@@ -12,6 +12,7 @@
#include "poly.h"
#include "io.h"
#include "sh.h"
#include "shcomplex.h"
ts_vesicle *initial_distribution_dipyramid(ts_uint nshell, ts_uint ncmax1, ts_uint ncmax2, ts_uint ncmax3, ts_double stepsize){
    ts_fprintf(stdout,"Starting initial_distribution on vesicle with %u shells!...\n",nshell);
@@ -94,7 +95,7 @@
    vesicle->pressure= tape->pressure;
    vesicle->pswitch=tape->pswitch;
    if(tape->shc>0){
        vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
        vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
    }
    else {
        vesicle->sphHarmonics=NULL;
src/main.c
@@ -12,6 +12,7 @@
#include "timestep.h"
#include "poly.h"
#include "sh.h"
#include "shcomplex.h"
/** Entrance function to the program
  * @param argv is a number of parameters used in program call (including the program name
@@ -53,7 +54,7 @@
        vesicle->clist->dmin_interspecies = tape->dmin_interspecies*tape->dmin_interspecies;
        /* spherical harmonics */
        if(tape->shc>0){
            vesicle->sphHarmonics=sph_init(vesicle->vlist,tape->shc);
            vesicle->sphHarmonics=complex_sph_init(vesicle->vlist,tape->shc);
        }
        else {
            vesicle->sphHarmonics=NULL;
src/shcomplex.c
New file
@@ -0,0 +1,127 @@
#include<math.h>
#include<stdlib.h>
#include<gsl/gsl_complex.h>
#include<gsl/gsl_complex_math.h>
#include<gsl/gsl_sf_legendre.h>
#include "general.h"
#include "sh.h"
#include "shcomplex.h"
ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l){
    ts_uint j,i;
    ts_spharm *sph=(ts_spharm *)malloc(sizeof(ts_spharm));
    sph->N=0;
    /* lets initialize Ylm for each vertex. */
    sph->Ylmi=(ts_double ***)calloc(l,sizeof(ts_double **));
    for(i=0;i<l;i++){
            sph->Ylmi[i]=(ts_double **)calloc(2*i+1,sizeof(ts_double *));
            for(j=0;j<(2*i+1);j++){
                sph->Ylmi[i][j]=(ts_double *)calloc(vlist->n,sizeof(ts_double));
            }
    }
    /* lets initialize ulm */
    sph->ulm=(ts_double **)calloc(l,sizeof(ts_double *));
    sph->ulmComplex=(gsl_complex **)calloc(l,sizeof(gsl_complex *));
    for(j=0;j<l;j++){
        sph->ulm[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
        sph->ulmComplex[j]=(gsl_complex *)calloc(2*j+1,sizeof(gsl_complex));
    }
    /* lets initialize sum of Ulm2 */
    sph->sumUlm2=(ts_double **)calloc(l,sizeof(ts_double *));
    for(j=0;j<l;j++){
        sph->sumUlm2[j]=(ts_double *)calloc(2*j+1,sizeof(ts_double));
    }
    /* lets initialize co */
//NOTE: C is has zero based indexing. Code is imported from fortran and to comply with original indexes we actually generate one index more. Also second dimension is 2*j+2 instead of 2*j+2. elements starting with 0 are useles and should be ignored!
    sph->co=(ts_double **)calloc(l+1,sizeof(ts_double *));
    for(j=0;j<=l;j++){
        sph->co[j]=(ts_double *)calloc(2*j+2,sizeof(ts_double));
    }
    sph->l=l;
    /* Calculate coefficients that will remain constant during all the simulation */
   precomputeShCoeff(sph);
    return sph;
}
ts_bool complex_sph_free(ts_spharm *sph){
    int i,j;
    if(sph==NULL) return TS_FAIL;
    for(i=0;i<sph->l;i++){
        if(sph->ulm[i]!=NULL) free(sph->ulm[i]);
        if(sph->ulmComplex[i]!=NULL) free(sph->ulmComplex[i]);
        if(sph->sumUlm2[i]!=NULL) free(sph->sumUlm2[i]);
        if(sph->co[i]!=NULL) free(sph->co[i]);
    }
        if(sph->co[sph->l]!=NULL) free(sph->co[sph->l]);
    if(sph->co != NULL) free(sph->co);
    if(sph->ulm !=NULL) free(sph->ulm);
    if(sph->ulmComplex !=NULL) free(sph->ulmComplex);
        if(sph->Ylmi!=NULL) {
            for(i=0;i<sph->l;i++){
                if(sph->Ylmi[i]!=NULL){
                    for(j=0;j<i*2+1;j++){
                        if(sph->Ylmi[i][j]!=NULL) free (sph->Ylmi[i][j]);
                    }
                    free(sph->Ylmi[i]);
                }
            }
            free(sph->Ylmi);
        }
    free(sph);
    return TS_SUCCESS;
}
ts_bool calculateUlmComplex(ts_vesicle *vesicle){
    ts_int i,j,k,m,l;
    ts_vertex *cvtx;
    ts_coord coord;
/* set all values to zero */
    for(i=0;i<vesicle->sphHarmonics->l;i++){
        for(j=0;j<2*i+1;j++) GSL_SET_COMPLEX(&(vesicle->sphHarmonics->ulmComplex[i][j]),0.0,0.0);
    }
    for(k=0;k<vesicle->vlist->n; k++){
        cvtx=vesicle->vlist->vtx[k];
    cart2sph(&coord,cvtx->x,cvtx->y,cvtx->z);
        for(i=0;i<vesicle->sphHarmonics->l;i++){
            for(j=0;j<2*i+1;j++){
        m=j-i;
        l=i;
        if(m>=0){
    //    fprintf(stderr, "Racunam za l=%d, m=%d\n", l,m);
                vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*gsl_sf_legendre_sphPlm(l,m,cos(coord.e3)))) );
        } else {
    //    fprintf(stderr, "Racunam za l=%d, abs(m=%d)\n", l,m);
                vesicle->sphHarmonics->ulmComplex[i][j]=gsl_complex_add(vesicle->sphHarmonics->ulmComplex[i][j], gsl_complex_conjugate(gsl_complex_mul_real(gsl_complex_polar(1.0,(ts_double)m*coord.e2),cvtx->solAngle*cvtx->relR*pow(-1,m)*gsl_sf_legendre_sphPlm(l,-m,cos(coord.e3)))) );
        }
            }
        }
    }
    return TS_SUCCESS;
}
ts_bool storeUlmComplex2(ts_vesicle *vesicle){
    ts_spharm *sph=vesicle->sphHarmonics;
    ts_int i,j;
    for(i=0;i<sph->l;i++){
            for(j=0;j<2*i+1;j++){
                sph->sumUlm2[i][j]+=gsl_complex_abs2(sph->ulmComplex[i][j]);
            }
    }
    sph->N++;
    return TS_SUCCESS;
}
src/shcomplex.h
New file
@@ -0,0 +1,8 @@
#ifndef _H_SH_COMPLEX
#define _H_SH_COMPLEX
#include "general.h"
ts_bool storeUlmComplex2(ts_vesicle *vesicle);
ts_spharm *complex_sph_init(ts_vertex_list *vlist, ts_uint l);
ts_bool complex_sph_free(ts_spharm *sph);
ts_bool calculateUlmComplex(ts_vesicle *vesicle);
#endif
src/tape
@@ -44,15 +44,15 @@
####### Program Control ############
#how many MC sweeps between subsequent records of states to disk
mcsweeps=10
mcsweeps=250
#how many initial mcsweeps*inititer MC sweeps before recording to disk?
inititer=0
#how many records do you want on the disk iteration are there in a run?
iterations=10
iterations=1000
###### Spherical harmonics ###########
spherical_harmonics_coefficients=21
spherical_harmonics_coefficients=15
#shut up if we are using cluster!!!
quiet=false
src/timestep.c
@@ -10,6 +10,7 @@
#include "io.h"
#include "stats.h"
#include "sh.h"
#include "shcomplex.h"
#include "vesicle.h"
ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations, ts_uint start_iteration){
@@ -51,9 +52,9 @@
            r0=getR0(vesicle);
            if(vesicle->sphHarmonics!=NULL){
                preparationSh(vesicle,r0);
                calculateYlmi(vesicle);
                calculateUlm(vesicle);
                storeUlm2(vesicle);
                //calculateYlmi(vesicle);
                calculateUlmComplex(vesicle);
                storeUlmComplex2(vesicle);
                saveAvgUlm2(vesicle);
                fd1=fopen("state.dat","w");
                fprintf(fd1,"%e %e\n",vesicle->volume, getR0(vesicle));
src/vesicle.c
@@ -7,6 +7,7 @@
#include "stdlib.h"
#include "poly.h"
#include "sh.h"
#include "shcomplex.h"
ts_vesicle *init_vesicle(ts_uint N, ts_uint ncmax1, ts_uint ncmax2, ts_uint
ncmax3, ts_double stepsize){
@@ -36,7 +37,7 @@
    triangle_list_free(vesicle->tlist);
    cell_list_free(vesicle->clist);
    poly_list_free(vesicle->poly_list);
    sph_free(vesicle->sphHarmonics);
    complex_sph_free(vesicle->sphHarmonics);
    free(vesicle);
    return TS_SUCCESS;
}