calc_vesicle_radius.cpp

#include <iostream>
#include <cassert>
#include <string>
#include <vector>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include "mvector3.hpp"

using namespace std;

const double cutof_leng = 1.5;
enum{
  CLUSTER_N = 2,
  MOL_NUMB  = 4,
};

double3 grid_leng(cutof_leng);
int grid_numb[3] = { -1 };
int all_grid_numb( -1 );
double3 L( -1.0 );

void LoadParticleDat(ifstream& fin, vector<double3>& pos, vector<float>& n_water, const int sys_size){
  float buf_f; double3 buf_d3;
  for(int i=0; i<sys_size; i++){
    fin >> buf_f >> buf_d3.x >> buf_d3.y >> buf_d3.z;
    n_water.push_back(buf_f); pos.push_back(buf_d3);
  }
}

void LoadScale(ifstream& fin){
  double buf[6] = {0.0};
  fin >> L.x >> L.y >> L.z >> buf[0] >> buf[1] >> buf[2] >> buf[3] >> buf[4] >> buf[5];
}

void SetScale(){
  grid_numb[0] = static_cast<int>(L.x / cutof_leng);
  grid_numb[1] = static_cast<int>(L.y / cutof_leng);
  grid_numb[2] = static_cast<int>(L.z / cutof_leng);
  all_grid_numb = grid_numb[0] * grid_numb[1] * grid_numb[2];
  grid_leng.x = L.x / grid_numb[0];
  grid_leng.y = L.y / grid_numb[1];
  grid_leng.z = L.z / grid_numb[2];
}

int GenHash(const double3& r){
  int i = (int)(r.x / grid_leng.x);
  int j = (int)(r.y / grid_leng.y);
  int k = (int)(r.z / grid_leng.z);
 
  if(i == grid_numb[0]) i--;
  if(j == grid_numb[1]) j--;
  if(k == grid_numb[2]) k--;
  
  const int hash = i + grid_numb[0] * (j + k * grid_numb[1]);
  assert(hash < all_grid_numb);
  assert(hash >= 0);
  return hash;
}

int GenHash(int* q){
  if(q[0] < 0 || q[0] > grid_numb[0]) q[0] -= (int)q[0]*grid_numb[0]/fabs(q[0]);
  if(q[1] < 0 || q[1] > grid_numb[1]) q[1] -= (int)q[1]*grid_numb[1]/fabs(q[1]);
  if(q[2] < 0 || q[2] > grid_numb[2]) q[2] -= (int)q[2]*grid_numb[2]/fabs(q[2]);
  const int hash = q[0] + grid_numb[0]*(q[1] + q[2]*grid_numb[1]);
  return hash;
}

void MinImage(double3& dr){
  dr.x -= L.x * round(dr.x / L.x);
  dr.y -= L.y * round(dr.y / L.y);
  dr.z -= L.z * round(dr.z / L.z);  
}

void GetNearImage(double& dr, double& rj, double L){
  if(dr < -0.5 * L){
    rj += L; dr += L;
  }else if(dr > 0.5 * L){
    rj -= L; dr -= L;
  }
}

void MinImage(double3& dr, double3& rj){
  GetNearImage(dr.x, rj.x, L.x);
  GetNearImage(dr.y, rj.y, L.y);
  GetNearImage(dr.z, rj.z, L.z);
}

void CreateNearlist(vector<vector<int> > &ilist){
  ilist.resize(all_grid_numb);
  for(int i=0; i<all_grid_numb; i++){
    ilist[i].resize(27);
  }
  
  int i_num=0;
  for(int iz=0; iz<grid_numb[2]; iz++){
    for(int iy=0; iy<grid_numb[1]; iy++){
      for(int ix=0; ix<grid_numb[0]; ix++){
	int j_num=0;
	for(int jz=-1; jz<2; jz++){
	  for(int jy=-1; jy<2; jy++){
	    for(int jx=-1; jx<2; jx++){
	      int q[3] = {ix+jx, iy+jy, iz+jz};
	      ilist[i_num][j_num]=GenHash(q);
	      j_num++;
	    }
	  }
	}
	i_num++;
      }
    }
  }
}

void RegistPrtclIdx(vector<int> &prtcl_idx, const vector<double3> &pos, vector<int> &elem_in_grid, vector<int> &adrs_grid, vector<int> &prtcl_in_cell){
  elem_in_grid.assign(elem_in_grid.size(), 0);
  adrs_grid[0] = 0;
  for(size_t i=0; i<pos.size(); i++){
    const int hash = GenHash(pos[i]);
    prtcl_in_cell[i] = hash;
    elem_in_grid[hash]++;
  }
  for(int i=0; i<all_grid_numb; i++){
    adrs_grid[i+1] = adrs_grid[i] + elem_in_grid[i];
  }
  for(size_t i=0; i<pos.size(); i++){
    const int hash = prtcl_in_cell[i];
    const int temp_idx = adrs_grid[hash];
    prtcl_idx[temp_idx] = i;
    adrs_grid[hash]++;
  }
  for(int i=0; i<all_grid_numb; i++){
    adrs_grid[i] = adrs_grid[i] - elem_in_grid[i];
  }
}

void CheckReg(const vector<int>& adrs_grid, const vector<int> &prtcl_idx, const vector<double3> &pr){
  for(int tar=0; tar<all_grid_numb; tar++){
    for(int pi = adrs_grid[tar]; pi<adrs_grid[tar+1]; pi++){
      const int hash = GenHash(pr[prtcl_idx[pi]]);
      assert(hash == tar);
    }
  }
}

void Dfs(int pid, vector<bool>& reg_flag, const vector<vector<int> >& n_prtcl_id, vector<int>& patch_id, vector<double3>& pos, int cur_patch){
  if(reg_flag[pid]) return;
  reg_flag[pid] = true;
  patch_id[pid] = cur_patch;
  const int n_prtcl_size = n_prtcl_id[pid].size();
  for(int i=0; i<n_prtcl_size; i++){
    const int	pjd = n_prtcl_id[pid][i];
    double3	dr  = pos[pjd] - pos[pid];
    MinImage(dr, pos[pjd]);
    Dfs(pjd, reg_flag, n_prtcl_id, patch_id, pos, cur_patch);
  }
}

int CountRegN(const vector<bool>& reg_flag){
  int n=0;
  for(size_t i=0; i<reg_flag.size(); i++) n += reg_flag[i];
  return n;
}

int SearchNextPid(const vector<bool>& reg_flag, size_t size){
  bool flag = true;
  int ret = 0;
  while(flag){
    int pid = rand() % size;
    if(!reg_flag[pid]){
      ret = pid;
      flag = false;
    }
  }
  return ret;
}

double3 CalcCMpos(const vector<double3>& pos){
  double3 cmpos(0.0);
  for(size_t i=0; i<pos.size(); i++){
    cmpos += pos[i];
  }
  cmpos /= pos.size();
  return cmpos;
}

void DumpRadiusDens(ofstream& fout0, ofstream& fout1, const vector<double3>& pos, const vector<float>& n_water, const vector<int> patch_id, const int cur_patch){
  if(cur_patch != 1) return;
  const double3 cm_pos = CalcCMpos(pos);
  struct rad_elem{
    double rad[2] = {0.0}; int elem[2] = {0};
  } r_elem;
  
  for(size_t i=0; i<pos.size(); i++){
    const double3 dr = pos[i] - cm_pos;
    const double rad = dr.norm2();
    const int   p_id = patch_id[i];
    r_elem.rad[p_id]  += rad;
    r_elem.elem[p_id] ++;
    fout0 << rad << " " << n_water[i] << endl;
  }

  r_elem.rad[0] /= r_elem.elem[0];
  r_elem.rad[1] /= r_elem.elem[1];

  if(r_elem.rad[0] < r_elem.rad[1]){
    fout1 << r_elem.rad[0] << " " << r_elem.rad[1] << endl;
  }else{
    fout1 << r_elem.rad[1] << " " << r_elem.rad[0] << endl;
  }
}

void DumpInOutElem(const vector<int>& patch_id, ofstream& fout, const int cur_patch, const int cur_time){
  if(cur_patch != CLUSTER_N) return;
  fout << cur_time << " ";
  int ret[CLUSTER_N] = {0};
  for(int i=0; i<cur_patch; i++) ret[i] = count(patch_id.begin(), patch_id.end(), i);
  sort(ret, ret + CLUSTER_N);
  for(int i=0; i<CLUSTER_N; i++) fout << ret[i] << " ";
  fout << endl;
}

void AssertFileOpen(const istream& ist, const string& str){
  if(!ist){
    cerr << "Cannot open " << str << endl;
    assert(ist);
  }
}

int main(int argc,char* argv[]){
  assert(argc == 3);
  srand((unsigned)time(NULL));
  const int beg_time = atoi(argv[2]);
  string s_cur  = argv[1];
  const string fname[5] = {
    s_cur + "/near_waters.txt",
    s_cur + "/macro_data.txt",
    s_cur + "/inout_elem.txt",
    s_cur + "/rad_water_dense.txt",
    s_cur + "/mean_rad.txt"
  };
  
  ifstream f_p(fname[0].c_str());
  ifstream f_m(fname[1].c_str());
  ofstream fout0(fname[2].c_str());
  ofstream fout1(fname[3].c_str());
  ofstream fout2(fname[4].c_str());
  
  int wN, hbN, all_time, time_step;
  double buf;

  AssertFileOpen(f_p, fname[0]);
  AssertFileOpen(f_m, fname[1]);
  
  //NOTE: sys_size = amphiphile molecular number
  f_m >> wN >> hbN >> L.x >> buf >> all_time >> time_step;
  const int sys_size = hbN / MOL_NUMB; 
  L.y = L.x;
  L.z = L.x;

  SetScale();

  //nearest cell
  vector<vector<int> > ilist;
  
  vector<int> elem_in_grid(2 * all_grid_numb, 0);
  vector<int> adrs_grid(2 * all_grid_numb + 1, 0);
  vector<double3> pos;
  vector<float>   n_water;
  vector<int    > prtcl_idx;
  vector<int    > prtcl_in_cell;
  vector<bool   > reg_flag;
  vector<vector<int> > n_prtcl_id;
  vector<int    > patch_id;
  
  for(int time=0; time < all_time; time += time_step){
    cout << time << endl;
    if(time >= beg_time){
      LoadParticleDat(f_p, pos, n_water, sys_size);
      CreateNearlist(ilist);

      const size_t size = pos.size();
      prtcl_in_cell.resize(size,-1);
      prtcl_idx.resize(size,-1);
      patch_id.resize(size,-1);
      n_prtcl_id.resize(size);
    
      RegistPrtclIdx(prtcl_idx, pos, elem_in_grid, adrs_grid, prtcl_in_cell);

      //debug
      CheckReg(adrs_grid, prtcl_idx, pos);

      //generate graph
      for(size_t idx=0; idx<size; idx++){
	const int pi     = prtcl_idx[idx];
	const int hash   = prtcl_in_cell[pi];
	const double3 ri = pos[pi];
	for(int i=0; i<27; i++){
	  const int intr_grid = ilist[hash][i];
	  const int beg_id    = adrs_grid[intr_grid    ];
	  const int end_id    = adrs_grid[intr_grid + 1];
	  for(int id=beg_id; id < end_id; id++){
	    const int pj = prtcl_idx[id];
	    double3 rj = pos[pj];
	    double3 dr = rj - ri;
	    MinImage(dr);
	    const double dist = dr.norm2();
	    if((dist < cutof_leng) && (pi != pj)){
	      n_prtcl_id[pi].push_back(pj);
	    }
	  }
	}
      }
      
      reg_flag.resize(size, false);
      int   cur_patch = 0;
      size_t regist_n = 0;
      int   pid       = 0;
      while(regist_n != size){
	pid = SearchNextPid(reg_flag, size);
	Dfs(pid, reg_flag, n_prtcl_id, patch_id, pos, cur_patch);
	regist_n = CountRegN(reg_flag);
	cout << regist_n << endl;
	cur_patch++;
      }

      DumpInOutElem(patch_id, fout0, cur_patch, time);
      DumpRadiusDens(fout1, fout2, pos, n_water, patch_id, cur_patch);
    }
    
    pos.clear();
    n_water.clear();
    prtcl_idx.clear();
    reg_flag.clear();
    for(size_t i=0; i<n_prtcl_id.size(); i++) n_prtcl_id[i].clear();
  }
  
  return 0;
}