タイトルのとおりです。ソースコード（md_water_ewald.cpp）は
https://github.com/poinco-gogo/md_water
で配布しています。学生時代に作ったものなので色々とツッコミどころ満載ですが、動くことは動きます（重要）。
力場は TIP3P で、 NVE/NVT のみの対応となります。
位置 Verlet 法を採用¹し、強制的に SHAKE します。²
温度制御は Langevin 方程式に従っています（BBK integrator）。³
式を explicit に書くなら

r(t+\Delta{t})=\frac{1}{B}\left\{2r(t)-Ar(t-\Delta{t})+\frac{\Delta t^2}{m}\left[F(t)+R(t)\right]\right\}\\
A=1-\gamma\frac{\Delta t}{2}\\
B=1+\gamma\frac{\Delta t}{2}\\
<R(0)R(t)>=2k_{B}T\gamma m \delta(t)

$R(t)$ はガウスランダム過程で、上式の最終行を満たします。 $T$ は熱浴の温度、 $\delta(t)$ はディラックのデルタ関数です。摩擦係数 $\gamma=0$ とすると通常の位置 Verlet の式に戻ることに注意してください。インプットの config ファイルで gamma 0とすることで NVE のシミュレーションをすることができます。

以下、コードの解説をしていきたいと思います。
（commit: 8b01725 に基づく）

冒頭

#include "vector3.h"

このプログラムでは 3 次元ベクトルの計算に
http://cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/SRC/vector3.h
を使用します（ただし、 abs を vabs に変更して下さい）。今なら Eigen を使う方があとで行列演算をしたくなった時に便利でしょう。

クラス定義（Atom）

class Atom
{       
        public: 
        Vector3 position;
        Vector3 velocity;
        Vector3 fold, fnew, vnew, rold, rnew;
        string  atomname;
        double  charge, mass, R; 
        Atom(double x, double y, double z)
        {       
                position.x = x;
                position.y = y;
                position.z = z;
        }
};

原子を記述するクラスです。R は温度制御の時に使用する定数です（下記参照）。

クラス定義（System）

class System
{
        public:
        int nstep, natom, kmax;
        double dt, T, gamma, boxsize, cutoff, volume, ewcoeff, lj, es;
        Vector3 origin, a1, a2, a3, g1, g2, g3, b1, b2, b3;

        System()
        {
                dt = 9999;
                T = 9999;
                gamma = 9999;
                nstep = 0;
                natom = 0;
                boxsize = 1;
                origin.x = 9999; origin.y = 9999; origin.z = 9999;
                cutoff = 1.;
                lj = 0;
                es = 0;
        }

        void setup_box()
        {
                volume = boxsize * boxsize * boxsize;
                a1.x = boxsize; a1.y = 0; a1.z = 0;
                a2.x = 0; a2.y = boxsize; a2.z = 0;
                a3.x = 0; a3.y = 0; a3.z = boxsize;
                double rvol2pi = 2. * M_PI / volume;
                g1 = (a2 % a3) * rvol2pi;
                g2 = (a3 % a1) * rvol2pi;
                g3 = (a1 % a2) * rvol2pi;
        }
};

系の情報を扱うクラスです。gamma は Langevin 方程式の friction coefficient というやつです。

main 関数

if (argc < 4)
{
        cout << "usage: ./a.out trjout psf pdb config\n";
        return 1;
}

引数が足りないときに使い方を表示します（return 1 でいいのか？）。

引数	意味
trjout	トラジェクトリを保存するアウトプットファイル名。トラジェクトリは XYZ 形式で保存します。
psf	系に存在するすべての分子の結合（bond, angle,...）の情報を記述するファイル（psf 形式）
pdb	初期座標（PDB 形式）
config	MD のパラメータを記述するファイル

Github に、 psf, pdb, config ファイルのサンプルがあります。

srand( (unsigned int)time(NULL) );

現在時刻をもとに乱数の発生系列を定義します。あとで rand() によって乱数を発生させていますが、 C++11 ではメルセンヌ・ツイスターが標準ライブラリに実装されているので、そちらを使う方がよいでしょう（参考）。

const int nfree = natom * 3 - nwat * 3;

系の自由度を定義しています。強制的に SHAKE しているので、水の結合の自由度を全体から引いておきます。

const double ps2asu = 20.455;

本プログラムでは原子の質量を原子質量単位（atomic mass unit）、エネルギーの単位を kcal/mol、長さの単位をオングストロームで表現しています。この場合、時間の単位はピコ秒に ps2asu をかけたものになります。

const double dt_ps = system.dt * 0.001; // in ps

インプットの config ファイルに dt 2 と書いた場合、それはタイムステップが 2 フェムト秒であることを意味します。なのでいったんピコ秒に変換します。

const double A = 1. - gamma * dt * 0.5;
const double B = 1. + gamma * dt * 0.5;
const double inB = 1. / B;

冒頭の式の $A$ と $B$ を計算しています。

for (int i = 0; i < natom; i++)
{      
        Atom& at = atomVector[i];
        at.R = sqrt(2. * kb * T * gamma * at.mass / dt);
}

各原子ごとに上式の $R$ の平均的な大きさを計算します。これに正規分布 $N(0, 1)$ に従う乱数を乗じたものが、ランダム力になります。

const double rOH  = 0.9572;
const double aHOH = 104.52 / 180. * acos(-1.0);

TIP3P のジオメトリを定義しています。

// make reciprocal vectors
vector<Vector3> g;
int kmax = system.kmax;
int sqkmax = kmax - 1;
sqkmax *= sqkmax;
double dum = 0;
for (int k = 0; k < kmax; k++)
for (int i = -kmax + 1; i < kmax; i++)
for (int j = -kmax + 1; j < kmax; j++)
{
        if (k * k + j * j + i * i > sqkmax or (i==0 && j==0 && k==0))
                continue;
        Vector3 vtmp(i * system.g1.x, j * system.g2.y, k * system.g3.z);
        if (vabs(vtmp) > dum)
                dum = vabs(vtmp);
        g.push_back(vtmp);
}

逆格子ベクトル $\vec{G}=2 \pi \left(\frac{k_{x}}{L_{x}},\frac{k_{y}}{L_{y}},\frac{k_{z}}{L_{z}} \right)$ を定義しています。 $k_{x}$, $k_{y}$ については、 config ファイルで指定した kmax の値をもとに -kmax から +kmax までとりますが、 $k_{z}$ については0 から +kmax までしかとりません。この条件の下で結果を 2 倍すれば求めるべき値を得ることができます。⁴
kmax は、静電相互作用エネルギーが収束するまで充分大きくとりましょう。
（エワルド係数 ewcoeff との兼ね合いで設定するべき kmax の値が変わることに注意してください。）

double ew_self = 0;
for (int i = 0; i < atomVector.size(); i++)
{
        ew_self += atomVector[i].charge * atomVector[i].charge;
}
ew_self *= -system.ewcoeff / sqrt(M_PI) * 332.0636;

エワルドの自己相互作用エネルギー項です。332.0636 は、$e^2/4 \pi \epsilon_{0}$ だったような気がしますが、記憶が定かではありません。

for (int i = 0; i < natom; i++)
{
        Atom& at = atomVector[i];
        at.velocity.x = gauss() * sqrt(kb * T / at.mass);
        at.velocity.y = gauss() * sqrt(kb * T / at.mass);
        at.velocity.z = gauss() * sqrt(kb * T / at.mass);
}

初速度をボルツマン分布からサンプリングしています。

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.fold = at.fnew;
        at.rold = at.position;
        at.position = at.rold + dt * at.velocity + dtdt_div2 / at.mass * at.fold;
}

冒頭の式で、 $r(t+\Delta t)$ を計算するために $r(t)$ と $r(t-\Delta t)$ が必要です。ここではインプットの座標が $r(t-\Delta t)$ だと思うことにして、 $r(t)$ をテイラー展開の 2 次までとることで求めます。

r(t) = r(0)+v(0)\Delta t + \frac{\Delta t^2}{2m}F(0)

for (int istep = 1; istep <= nstep; istep++)
{

この中が MD のメインループです。

calc_frc(atomVector, lj_pair_list, el_pair_list, system, g);

まず position に保存された座標に基づいて force を計算します。計算された force は fnew に保存されます。

for (int i = 0; i < atomVector.size(); i++)
{
        Vector3 noise( gauss(), gauss(), gauss());
        Atom& at = atomVector[i];
        at.rnew = 2. * at.position - at.rold + gamma * dt_div2 * at.rold + dt * dt / at.mass * (at.fnew + at.R * noise);
        at.rnew = at.rnew * inB;
}

冒頭の式に基づいて $r(t+\Delta t)$ を計算します（ rnew に保存される）。

for (int ishake = 0; ishake < 100; ishake++)
{
        if (shake(atomVector, shake_list))
                break;
        if (ishake == 99)
        {
                cerr << "error: shake does not converged.\n";
                return 1;
        }
}

SHAKE ループです。

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.vnew = 0.5 / dt * (at.rnew - at.rold);
}

$v(t+\Delta t)=\frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}$ で、更新された速度ベクトルを計算します。結果は vnew に保存されます。

double K = calc_kin(atomVector);
calc_pot(atomVector, lj_pair_list, el_pair_list, system, g);
system.es += ew_self;

運動エネルギー、ポテンシャルエネルギーを計算します。
運動エネルギーは vnew に基づいて計算しています。
ポテンシャルエネルギーは position に基づいて計算しているみたいです。これは変ですね。 rnew に基づく計算であるべきではないだろうか？

for (int i = 0; i < atomVector.size(); i++)
{
        Atom& at = atomVector[i];
        at.rold = at.position;
        at.position = at.rnew;
        at.velocity = at.vnew;
}

次のステップに備えて変数の値を入れ替えています。

以上で MD のループは終わりです。
ちょっと長くなってきたので、続きは別にします。
続き→水（TIP3P）の古典 MD プログラム part 2

拘束の導入が簡単なので学習目的にはよいでしょう。しかし速度 Verlet の方が velocity のエラーが少ないので良いかもしれません。有名な MD プログラムは大抵速度 Verlet か leap frog なのではないでしょうか。 ↩
時間の刻み幅を大きくしたいからですね。 ↩
BBK integrator の速度 velocity バージョンは Bussi, G.; Parrinello, M. "Accurate sampling using Langevin dynamics", Phys. Rev. E, 2007, 75, 056707. で提案された integrator の、時間の刻み幅が小さい極限で一致する、気がする。 ↩
Tuckerman, M. E. "Statistical Mechanics: Theory and Molecular Simulation", Oxford Graduate Texts, pp. 661 ↩

水（TIP3P）の古典 MD プログラム part 1

冒頭

クラス定義（Atom）

クラス定義（System）

main 関数