Skip to content

Commit

Permalink
flat membrane code
Browse files Browse the repository at this point in the history
  • Loading branch information
vikash committed Jan 20, 2023
1 parent d079017 commit 1ab3b9d
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 36 deletions.
68 changes: 50 additions & 18 deletions mains/flat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,34 @@ void check_absurd_neighbours(MESH mesh, int N){

}
}

void debug_ener(Vec3d *pos, MESH mesh,
double *lij_t0, MBRANE_para mbrane,
MCpara mcpara, AFM_para afm,
ActivePara activity, SPRING_para spring){

int idx;
int cm_idx, num_nbr;
double E_b;

for(idx = 0; idx < mbrane.N; idx++){
cm_idx = mesh.nghst*idx;
num_nbr = mesh.numnbr[idx];


/* E_b = bending_energy_ipart_neighbour(pos, mesh, idx, mbrane); */
E_b = stretch_energy_ipart(pos, (int *) (mesh.node_nbr_list + cm_idx),
(double *) (lij_t0 + cm_idx), num_nbr,
idx, mbrane);


/* E_b = bending_energy_ipart(pos, */
/* (int *) (mesh.node_nbr_list + cm_idx), */
/* num_nbr, idx, mbrane); */

printf(" curvature %d %lf \n",idx, E_b);
}
}
//
int main(int argc, char *argv[]){
pid_t pid = getpid();
Expand All @@ -29,6 +57,7 @@ int main(int argc, char *argv[]){
Vec3d afm_force,spring_force[2];
FILE *fid;
double *lij_t0;
bool *is_attractive;
double Pole_zcoord;
string outfolder,syscmds, para_file, log_file, outfile, filename;
int mpi_err,mpi_rank=0;
Expand All @@ -43,7 +72,7 @@ int main(int argc, char *argv[]){
seed_v = 12343234;
init_rng(seed_v);
//
outfolder = ZeroPadNumber(mpi_rank)+"/";
outfolder = argv[1]; //Z"e"roPadNumber(mpi_rank)+"/";
cout << "I am in folder "+ outfolder << endl;
filename = outfolder + "/para_file.in";
write_param(outfolder + "/para.out",mbrane,mcpara,spring);
Expand All @@ -65,6 +94,7 @@ int main(int argc, char *argv[]){
/* define all the paras */
mbrane.tot_energy = (double *)calloc(1, sizeof(double));
mbrane.len = 1.0;
mbrane.is_fluid = true;
activity.activity = (double *)calloc(mbrane.N, sizeof(double));
mbrane.tot_energy[0] = 0e0;
init_activity(activity, mbrane.N);
Expand All @@ -81,16 +111,29 @@ int main(int argc, char *argv[]){
hdf5_io_read_mesh((int *) mesh.numnbr,
(int *) mesh.node_nbr_list, outfolder+"/nbrs.h5");
init_eval_lij_t0(Pos, mesh, lij_t0, &mbrane, &spring);

// Et[0] = stretch_energy_total(Pos, mesh, lij_t0, mbrane);



for(iter=1; iter <= mcpara.tot_mc_iter; iter++){
num_moves = monte_carlo_fluid(Pos, mesh,
mbrane, mcpara, afm, activity, spring);
check_absurd_neighbours(mesh, mbrane.N);

printf("%d %d \n", iter, num_moves );
num_moves = monte_carlo_3d(Pos, mesh, lij_t0, is_attractive,
mbrane, mcpara, afm, activity, spring);

printf("elastic part stats %d %d \n", iter, num_moves );
/* num_moves = monte_carlo_fluid(Pos, mesh, */
/* mbrane, mcpara, afm, activity, spring); */
/* printf("fluid stats %d %d \n", iter, num_moves ); */
if(i%mcpara.dump_skip == 0){
outfile=outfolder+"/snap_"+ZeroPadNumber(iter/mcpara.dump_skip)+".h5";
// sprintf(outfile,"%s/snap_%04d.h5",outfolder,(int)(i/mcpara.dump_skip));
hdf5_io_write_pos((double*) Pos, 3*mbrane.N, outfile);
/* syscmds="cp "+outfile+" "+outfolder+"/restart.h5"; */
system(syscmds.c_str());
}

/* check_absurd_neighbours(mesh, mbrane.N); */

}

// /************************************/
Expand All @@ -103,18 +146,7 @@ int main(int argc, char *argv[]){
// Et[0] = stretch_energy_total(Pos, mesh, lij_t0, mbrane);
// cout << "iter = " << i << "; Accepted Moves = " << (double) num_moves*100/mcpara.one_mc_iter << " %;"<<
// " totalener = "<< mbrane.tot_energy[0] << "; volume = " << vol_sph << endl;
// if(i%mcpara.dump_skip == 0){
// outfile=outfolder+"/snap_"+ZeroPadNumber(i/mcpara.dump_skip)+".h5";
// // sprintf(outfile,"%s/snap_%04d.h5",outfolder,(int)(i/mcpara.dump_skip));
// hdf5_io_write_pos((double*) Pos, 3*mbrane.N, outfile);
// syscmds="cp "+outfile+" "+outfolder+"/restart.h5";
// system(syscmds.c_str());
// }

// // num_moves = monte_carlo_3d(Pos, mesh, lij_t0, is_attractive,
// // mbrane, mcpara, afm, activity, spring);
// //
// num_moves = monte_carlo_fluid(Pos, mesh, mbrane, mcpara, afm, activity, spring);
// num_moves = monte_carlo_fluid(Pos, mesh, mbrane, mcpara, afm, activity, spring);

// }
/* fclose(fid); */
Expand Down
14 changes: 8 additions & 6 deletions src/Metropolis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -318,16 +318,18 @@ int monte_carlo_3d(Vec3d *pos, MESH mesh,
// de_vol = (2*dvol/(ini_vol*ini_vol))*(mbrane.volume[0] - ini_vol)
// + (dvol/ini_vol)*(dvol/ini_vol);
// de_vol = KAPPA*de_vol;

de = (Efin - Eini); //+ de_vol + de_pressure;
if (mcpara.algo == "mpolis"){
yes=Metropolis(de, activity.activity[idx], mcpara);
}else if(mcpara.algo == "glauber"){
yes=Glauber(de, activity.activity[idx], mcpara);
}

if (yes){
move = move + 1;
mbrane.tot_energy[0] += de;
mbrane.volume[0] += dvol;
/* mbrane.volume[0] += dvol; */
}else{
pos[idx].x = x_o;
pos[idx].y = y_o;
Expand All @@ -337,8 +339,7 @@ int monte_carlo_3d(Vec3d *pos, MESH mesh,
return move;
}

int monte_carlo_surf2d(Vec2d *Pos,
Nbh_list *neib, LJpara para,
int monte_carlo_surf2d(Vec2d *Pos, Nbh_list *neib, LJpara para,
MCpara mcpara, char *metric){

/// @brief Monte-Carlo routine for the points on the surface of sphere or flat
Expand All @@ -350,6 +351,7 @@ int monte_carlo_surf2d(Vec2d *Pos,
/// @param metric Topology of the surface "cart" for flat plane "sph" for
/// sphere
/// @return number of accepted moves

int i, move;
double x_o, y_o, x_n, y_n;
double de, Eini, Efin;
Expand Down Expand Up @@ -519,11 +521,11 @@ int monte_carlo_fluid(Vec3d *pos, MESH mesh,
/* bef_ij = diff_pbc(pos[idx_del1] , pos[idx_del2], mbrane.len); */
/* aft_ij = diff_pbc(pos[idx_add1] , pos[idx_add2], mbrane.len); */

double de = (norm(bef_ij) - mbrane.av_bond_len)*(norm(bef_ij) - mbrane.av_bond_len) -
(norm(aft_ij) - mbrane.av_bond_len)*(norm(aft_ij) - mbrane.av_bond_len);
double de = (norm(aft_ij) - mbrane.av_bond_len)*(norm(aft_ij) - mbrane.av_bond_len) -
(norm(bef_ij) - mbrane.av_bond_len)*(norm(bef_ij) - mbrane.av_bond_len);

// de = (Efin - Eini) + de_vol + de_pressure;
de = -1e4*de;
de = mbrane.YY*de;
if (mcpara.algo == "mpolis"){
yes=Metropolis(de, 0.0e0, mcpara);
}else if(mcpara.algo == "glauber"){
Expand Down
12 changes: 3 additions & 9 deletions src/forces_surf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,13 +245,7 @@ double bending_energy_ipart(Vec3d *pos, int *node_nbr, int num_nbr,
xikp = pos[idx]- pos[kpdx];
xjkp = pos[jdx]- pos[kpdx];
//
// xij = Vec3d_add(pos[idx], pos[jdx], -1e0);
// xijp1 = Vec3d_add(pos[idx], pos[jdxp1], -1e0);
// xik = Vec3d_add(pos[idx], pos[kdx], -1e0);
// xjk = Vec3d_add(pos[jdx], pos[kdx], -1e0);
// xikp = Vec3d_add(pos[idx], pos[kpdx], -1e0);
// xjkp = Vec3d_add(pos[jdx], pos[kpdx], -1e0);
//
//
lijsq = inner_product(xij,xij);
liksq = inner_product(xik,xik);
ljksq = inner_product(xjk,xjk);
Expand Down Expand Up @@ -302,7 +296,7 @@ double bending_energy_ipart_neighbour(Vec3d *pos,

for (j = idx*mesh.nghst; j < idx*mesh.nghst + mesh.numnbr[idx]; j++){
nbr = mesh.node_nbr_list[j];
num_nbr_j = mesh.numnbr[j];
num_nbr_j = mesh.numnbr[nbr];
cm_idx_nbr = nbr*mesh.nghst;
be += bending_energy_ipart(pos,
(int *) mesh.node_nbr_list + cm_idx_nbr,
Expand All @@ -329,7 +323,7 @@ double bending_energy_ipart_neighbour(Vec3d *pos,

vol = 0e0;
for(idx = 0; idx < para.N; idx++){
// /* idx = 2; */
/* idx = 2; */
cm_idx = idx*mesh.nghst;
num_nbr = mesh.numnbr[idx];

Expand Down
8 changes: 5 additions & 3 deletions src/init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,8 +154,8 @@ void init_eval_lij_t0(Vec3d *Pos, MESH mesh, double *lij_t0,
cm_idx = mesh.nghst * i;
for(k = cm_idx; k < cm_idx + num_nbr; k++) {
j = mesh.node_nbr_list[k];
dr = diff_pbc(Pos[i] , Pos[j], para->len);
/* dr = Pos[j] - Pos[i]; */
/* dr = diff_pbc(Pos[i] , Pos[j], para->len); */
dr = Pos[j] - Pos[i];
lij_t0[k] = sqrt(dr.x*dr.x + dr.y*dr.y + dr.z*dr.z);
sum_lij += sqrt(dr.x*dr.x + dr.y*dr.y + dr.z*dr.z);
npairs++;
Expand All @@ -167,7 +167,9 @@ void init_eval_lij_t0(Vec3d *Pos, MESH mesh, double *lij_t0,
r0=para->av_bond_len;
spring->constant=para->coef_bend/(r0*r0);
if(para->is_fluid){
for(i = 0; i < mesh.nghst*para->N; i++) lij_t0[i] = para->av_bond_len;
for(i = 0; i < mesh.nghst*para->N; i++){
lij_t0[i] = para->av_bond_len;
}
}
}

Expand Down

0 comments on commit 1ab3b9d

Please sign in to comment.