Use of Ewald Summation in MD++

MD++: Use of the Ewald Summation

Eunseok Lee

Compile

In MD++, the BMB potential (short range interaction) + the Ewald summation (long range interaction) is implemented. You can download the source code from the repository of MD++, then use 'make' to compile the binary.

> make bmb

When you complie it please notice that the link to FFTw package is required.

Parameters

Following parameters are used.

MD++ Ewald_Rcut= 6              # set R_cut 
MD++ Ewald_CE_or_PME = 0        # 0:CE, 1:PME
MD++ Ewald_option_Alpha = 1     # Automatic generation of Alpha
#MD++ Ewald_volume_change = 1    # Consider the volume change; only used in PME
MD++ Ewald_precision = 4.05     # Ewald_Alpha = Ewald_precision/Ewald_Rcut
#MD++ Ewald_Alpha = 0.45         # if Ewald_option_Alpha=0, manual input of Ewald_Alpha is necessary
MD++ Ewald_time_ratio = 2.5
#MD++ PME_K1 =  75
#MD++ PME_K2 =  75
#MD++ PME_K3 =  75
#MD++ PME_bsp_n = 10
MD++ Ewald_init                 # Initialize Ewald

Check the accuracy

In principle, the Ewald summation must be same regardless of alpha, which just changes the speed of calculation. To check if the Ewald summation is being calculated right, we can print out the energy and the Ewald summation for different Rcut (or alpha). In the source code, ewald.cpp, recover the commented out part (line #424 ~ #436), and then use 'MD++ eval' in script for different alphas. Use script file, ysz_Rcut_test.tcl (copy this file as scripts/Work/ysz_Rcut_test.tcl), make a directory of 'runs/YSZ_Rcut_Test/', and then run bash script, autorun_ysz_Rcut_test (copy this file as autorun_ysz_Rcut_test and make it executable). The result is like following.

Rcut   EPOT       EPOT_BMB    EPOT_Ewald  Ewald_Real   Ewald_Rec   Ewald_Self
3     -11746.345  2211.16828 -13957.514   -0.9369493   14274.1834 -28230.760
3.5   -11746.901  2210.61300 -13957.514   -11.601195   10251.8818 -24197.794
4     -11757.810  2199.70358 -13957.513   -60.336868   7275.89362 -21173.070
4.5   -11749.271  2208.24295 -13957.514   -188.77258   5051.76529 -18820.507
5     -11735.388  2222.12548 -13957.513   -429.41484   3410.35737 -16938.456
5.5   -11736.027  2221.48641 -13957.513   -790.66055   2231.74371 -15398.596
6     -11736.624  2220.88943 -13957.513   -1256.9764   1414.84287 -14115.380
6.5   -11737.860  2219.65295 -13957.513   -1798.0317   870.099924 -13029.581
7     -11738.553  2218.96011 -13957.513   -2379.3621   520.745745 -12098.897
7.5   -11738.648  2218.86439 -13957.513   -2970.2139   305.004859 -11292.304
8     -11738.994  2218.51979 -13957.514   -3547.4128   176.433772 -10586.535
8.5   -11739.253  2218.26181 -13957.515   -4095.9922   102.274376 -9963.7978
9     -11739.382  2218.14006 -13957.522   -4607.9822   60.7131600 -9410.2535
9.5   -11739.487  2218.06221 -13957.549   -5080.5281   37.9553160 -8914.9770
10    -11739.670  2217.96166 -13957.632   -5514.0766   25.6724522 -8469.2282
10.5  -11739.910  2217.93692 -13957.847   -5910.9586   19.0430761 -8065.9316
11    -11740.415  2217.91699 -13958.332   -6274.4151   15.3808837 -7699.2983
11.5  -11741.411  2217.90135 -13959.312   -6608.0074   13.2410181 -7364.5462
12    -11743.213  2217.89723 -13961.110   -6915.2914   11.8712117 -7057.6901
12.5  -11746.260  2217.89241 -13964.152   -7199.6613   10.8912866 -6775.3825
13    -11751.070  2217.89136 -13968.962   -7464.2854   10.1140953 -6514.7909
13.5  -11758.246  2217.89125 -13976.138   -7712.0839   9.44826301 -6273.5023
14    -11768.429  2217.89085 -13986.320   -7945.7252   8.85377107 -6049.4487
14.5  -11782.276  2217.89084 -14000.167   -8167.6273   8.30735539 -5840.8470
15    -11800.421  2217.89084 -14018.312   -8379.9592   7.79897806 -5646.1521
15.5  -11823.447  2217.89084 -14041.338   -8584.6430   7.32242251 -5464.0181

From this result (EPOT_Ewald) we can know that the safe range of R_cut is under about 8.5(A). In this example the size of supercell is 3x3x3. It is known that R_cut=9(A) is enough for 6x6x6 supercell (Courtesy to Hark Lee).