title: How to calculate Mean-squared displacement from a file in C++?
body: <p>I want to calculate the mean-squared displacements of atoms from a file in C++. I have the file where the atoms positions in each time step is written like below where 1..5 column respectively represents timestep, atom, and x,y,z co-ordinates of that atom.</p>
  
  <pre><code>Timestep = 0000000
  
  0000000   Fe      0     0     0
  0000000   Fe      0     1     0
  ......    ..      ..    ..    ..
  
  Timestep = 1000000
  
  1000000   Fe      0     0     1
  1000000   Fe      0     0     2
  .......   ..      ..    ..    ..
  And so on...
  </code></pre>
  
  <p>What I want is the length of the vector difference of the atoms and the time difference and then calculate MSD in a separate file (say 'msd.txt') where the 1st column would be in this particular case 1000000 and the 2nd column would be ([001 - 000]^2 + [002 - 010]^2)/2 i.e., squared of the displacement of the atoms. The displacement of individual atom is simply the length of the vector difference of individual atom. As in this case the 1st atom moves to [001] position from [000] position. A simple C++ code would be really helpful in this regard as I am really new in C++. I have written a sample code for this as shown below: </p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;vector&gt;
  #include &lt;cmath&gt;
  
  using namespace std;
  
  int Natoms = 16; //Number of atoms
  
  
  double length() const  // computes the length of the vector
    {
      return sqrt(x * x + y * y + z * z);
    }
  
  
  int main ()
  {
    ofstream oFile("msd.txt");
    oFile &lt;&lt; "# time  msd" &lt;&lt; endl;
    vector&lt;vector&lt;double&gt; &gt; table;  // stores 'dr' for all the atoms and different time steps
     for (int t = 0; t &lt; MaxSteps; t++)
     {
          timestep = 1000
          ct = (t+1)*timestep;
          vector&lt;double&gt; column;
          vector&lt;double&gt; rt; vector&lt;double&gt; r0;
          for (int i = 0; i &lt; Natoms; i++)
          {
                  double dr = rt[iTarget] - r0[iOrigin];
                  double dr2 = dr * dr;
                  column.push_back(dr2);
          }
          double sum_of_dr2 = 0. ; double MSD = 0. ;
          for(std::vector&lt;double&gt;::iterator it = column.begin(); it != column.end(); ++it)
                  sum_of_dr2 += *it;
  
          MSD = sum_of_dr2 * (1./(double)Natoms);
          table.push_back(column);
          oFile &lt;&lt; ct &lt;&lt; "\t" &lt;&lt; MSD &lt;&lt; endl;
     }
    oFile.close();
  
    return 0;
  }
  </code></pre>
  
  <p>But I don't know how to read from the given file so that I can have the <code>length</code> of the vector difference for each atom stored in <code>dr</code> in the sample code. Here <code>length</code> of the vector of 1st and 2nd atom after a time step 1000000 would be sqrt((0-0)^2+(0-0)^2+(0-1)^2) and sqrt((0-0)^2+(0-0)^2+(0-2)^2) resp. These data should be stored in <code>dr[i]</code> (say) where <code>i</code> = 1,2,..N is the no of atom.</p>