C-Example1

C语言实例1:铝原子体系势能与体积的变化关系

申明:下文出现的所有C语言代码均来自网上,也不涉及版权问题。这里仅供学习之用。源码为纯C语言编写。

文件结构

1
2
3
4
5
initfcc.h
initfcc.c
alpotential.h
alpotential.c
MD_main.c

文件内容

  • initfcc.h文件,其中声明了铝原子的初始FCC晶体参数函数。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    /*
    initfcc.h

    Created by Anders Lindman on 2013-03-15.
    */

    #ifndef _initfcc_h
    #define _initfcc_h

    extern void init_fcc(double[][3], int, double);


    #endif
  • initfcc.c文件,定义了初始FCC晶体参数。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    /*
    initfcc.c
    Program that arranges atoms on a fcc lattice.
    Created by Anders Lindman on 2013-03-15.
    */

    #include <stdio.h>

    /* Function takes a matrix of size [4*N*N*N][3] as input and stores a fcc lattice in it. N is the number of unit cells in each dimension and lattice_param is the lattice parameter. */
    void init_fcc(double positions[][3], int N, double lattice_param)
    {
    int i, j, k;
    int xor_value;

    for (i = 0; i < 2 * N; i++){
    for (j = 0; j < 2 * N; j++){
    for (k = 0; k < N; k++){
    if (j % 2 == i % 2 ){
    xor_value = 0;
    }
    else {
    xor_value = 1;
    }
    positions[i * N * 2 * N + j * N + k][0] = lattice_param * (0.5 * xor_value + k);
    positions[i * N * 2 * N + j * N + k][1] = lattice_param * (j * 0.5);
    positions[i * N * 2 * N + j * N + k][2] = lattice_param * (i * 0.5);
    }
    }
    }
    }
  • alpotential.h文件,声明了得到铝原子的力,势能和维里函数。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    /*
    alpotential.h

    Created by Anders Lindman on 2013-03-15.
    */

    #ifndef _alpotential_h
    #define _alpotential_h

    extern void get_forces_AL(double[][3] , double[][3], double, int);
    extern double get_energy_AL(double[][3], double, int);
    extern double get_virial_AL(double[][3], double, int);


    #endif
  • alpotential.c文件,定义了得到铝原子体系中的力,势能和维里的函数。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    179
    180
    181
    182
    183
    184
    185
    186
    187
    188
    189
    190
    191
    192
    193
    194
    195
    196
    197
    198
    199
    200
    201
    202
    203
    204
    205
    206
    207
    208
    209
    210
    211
    212
    213
    214
    215
    216
    217
    218
    219
    220
    221
    222
    223
    224
    225
    226
    227
    228
    229
    230
    231
    232
    233
    234
    235
    236
    237
    238
    239
    240
    241
    242
    243
    244
    245
    246
    247
    248
    249
    250
    251
    252
    253
    254
    255
    256
    257
    258
    259
    260
    261
    262
    263
    264
    265
    266
    267
    268
    269
    270
    271
    272
    273
    274
    275
    276
    277
    278
    279
    280
    281
    282
    283
    284
    285
    286
    287
    288
    289
    290
    291
    292
    293
    294
    295
    296
    297
    298
    299
    300
    301
    302
    303
    304
    305
    306
    307
    308
    309
    310
    311
    312
    313
    314
    315
    316
    317
    318
    319
    320
    321
    322
    323
    324
    325
    326
    327
    328
    329
    330
    331
    332
    333
    334
    335
    336
    337
    338
    339
    340
    341
    342
    343
    344
    345
    346
    347
    348
    349
    350
    351
    352
    353
    354
    355
    356
    357
    358
    359
    360
    361
    362
    363
    364
    365
    366
    367
    368
    369
    370
    371
    372
    373
    374
    375
    376
    377
    378
    379
    380
    381
    382
    383
    384
    385
    386
    387
    388
    389
    390
    391
    392
    393
    394
    395
    396
    397
    398
    399
    400
    401
    402
    403
    404
    405
    406
    407
    408
    409
    410
    411
    412
    413
    414
    415
    416
    417
    418
    419
    420
    421
    422
    423
    424
    425
    426
    427
    428
    429
    430
    431
    432
    433
    434
    435
    436
    437
    438
    439
    440
    441
    442
    443
    444
    445
    446
    447
    448
    449
    450
    451
    452
    453
    454
    455
    456
    457
    458
    459
    460
    461
    462
    463
    464
    465
    466
    467
    468
    469
    470
    471
    472
    473
    /*
    alpotential.c
    Program that contains functions that calculate properties (potential energy, forces, etc.) of a set of Aluminum atoms using an embedded atom model (EAM) potential.
    Created by Anders Lindman on 2013-03-14.
    */

    //使用嵌入原子势计算铝原子体系中的原子势能和力 的函数

    //头文件
    #include <stdio.h> //输入输出
    #include <math.h> //数学库
    #include <stdlib.h> //标准库

    /*Parameters for the AL EAM potential */
    #define PAIR_POTENTIAL_ROWS 18 //对势能列

    //对势
    const double pair_potential[90] = {2.0210, 2.2730, 2.4953, 2.7177, 2.9400, 3.1623, 3.3847, 3.6070, 3.8293, 4.0517, 4.2740, 4.4963, 4.7187, 4.9410, 5.1633, 5.3857, 5.6080, 6.0630, 2.0051, 0.7093, 0.2127, 0.0202, -0.0386, -0.0492, -0.0424, -0.0367, -0.0399, -0.0574, -0.0687, -0.0624, -0.0492, -0.0311, -0.0153, -0.0024, -0.0002, 0, -7.2241, -3.3383, -1.3713, -0.4753, -0.1171, 0.0069, 0.0374, 0.0122, -0.0524, -0.0818, -0.0090, 0.0499, 0.0735, 0.0788, 0.0686, 0.0339, -0.0012, 0, 9.3666, 6.0533, 2.7940, 1.2357, 0.3757, 0.1818, -0.0445, -0.0690, -0.2217, 0.0895, 0.2381, 0.0266, 0.0797, -0.0557, 0.0097, -0.1660, 0.0083, 0, -4.3827, -4.8865, -2.3363, -1.2893, -0.2907, -0.3393, -0.0367, -0.2290, 0.4667, 0.2227, -0.3170, 0.0796, -0.2031, 0.0980, -0.2634, 0.2612, -0.0102, 0};


    #define ELECTRON_DENSITY_ROWS 15 //电子密度列

    //电子密度
    const double electron_density[75] = {2.0210, 2.2730, 2.5055, 2.7380, 2.9705, 3.2030, 3.4355, 3.6680, 3.9005, 4.1330, 4.3655, 4.5980, 4.8305, 5.0630, 6.0630, 0.0824, 0.0918, 0.0883, 0.0775, 0.0647, 0.0512, 0.0392, 0.0291, 0.0186, 0.0082, 0.0044, 0.0034, 0.0027, 0.0025, 0.0000, 0.0707, 0.0071, -0.0344, -0.0533, -0.0578, -0.0560, -0.0465, -0.0428, -0.0486, -0.0318, -0.0069, -0.0035, -0.0016, -0.0008, 0, -0.1471, -0.1053, -0.0732, -0.0081, -0.0112, 0.0189, 0.0217, -0.0056, -0.0194, 0.0917, 0.0157, -0.0012, 0.0093, -0.0059, 0, 0.0554, 0.0460, 0.0932, -0.0044, 0.0432, 0.0040, -0.0392, -0.0198, 0.1593, -0.1089, -0.0242, 0.0150, -0.0218, 0.0042, 0};

    #define EMBEDDING_ENERGY_ROWS 13 //嵌入能列

    //嵌入能
    const double embedding_energy[65] = {0, 0.1000, 0.2000, 0.3000, 0.4000, 0.5000, 0.6000, 0.7000, 0.8000, 0.9000, 1.0000, 1.1000, 1.2000, 0, -1.1199, -1.4075, -1.7100, -1.9871, -2.2318, -2.4038, -2.5538, -2.6224, -2.6570, -2.6696, -2.6589, -2.6358, -18.4387, -5.3706, -2.3045, -3.1161, -2.6175, -2.0666, -1.6167, -1.1280, -0.4304, -0.2464, -0.0001, 0.1898, 0.2557, 86.5178, 44.1632, -13.5018, 5.3853, -0.3996, 5.9090, -1.4103, 6.2976, 0.6785, 1.1611, 1.3022, 0.5971, 0.0612, -141.1819, -192.2166, 62.9570, -19.2831, 21.0288, -24.3978, 25.6930, -18.7304, 1.6087, 0.4704, -2.3503, -1.7862, -1.7862};


    /* Evaluates the spline in x. */

    double splineEval(double x, const double *table,int m) {
    /* int m = mxGetM(spline), i, k;*/
    int i, k;

    /*double *table = mxGetPr(spline);*/
    double result;

    int k_lo = 0, k_hi = m;

    /* Find the index by bisection. */
    while (k_hi - k_lo > 1) {
    k = (k_hi + k_lo) >> 1; //右移运算,二进制右移1,即偶数数值变为原来的一半,奇数减一除以二。
    if (table[k] > x)
    k_hi = k;
    else
    k_lo = k;
    }

    /* Switch to local coord. */
    x -= table[k_lo];

    /* Horner's scheme */
    result = table[k_lo + 4*m];
    for (i = 3; i > 0; i--) {
    result *= x;
    result += table[k_lo + i*m];
    }

    return result;
    }

    /* Evaluates the derivative of the spline in x. */

    double splineEvalDiff(double x, const double *table, int m) {
    /*int m = mxGetM(spline), i, k;
    double *table = mxGetPr(spline);
    */
    int i, k;
    double result;

    int k_lo = 0, k_hi = m;

    /* Find the index by bisection. */
    while (k_hi - k_lo > 1) {
    k = (k_hi + k_lo) >> 1;
    if (table[k] > x)
    k_hi = k;
    else
    k_lo = k;
    }

    /* Switch to local coord. */
    x -= table[k_lo];

    /* Horner's scheme */
    result = 3*table[k_lo + 4*m];
    for (i = 3; i > 1; i--) {
    result *= x;
    result += (i-1)*table[k_lo + i*m];
    }

    return result;
    }

    /* Returns the forces */
    //返回力
    void get_forces_AL(double forces[][3], double positions[][3], double cell_length, int nbr_atoms)
    {
    int i, j;
    double cell_length_inv, cell_length_sq;
    double rcut, rcut_sq;
    double densityi, dens, drho_dr, force;
    double dUpair_dr;
    double sxi, syi, szi, sxij, syij, szij, rij, rij_sq;
    //分配内存空间
    double *sx = malloc(nbr_atoms * sizeof (double));
    double *sy = malloc(nbr_atoms * sizeof (double));
    double *sz = malloc(nbr_atoms * sizeof (double));
    double *fx = malloc(nbr_atoms * sizeof (double));
    double *fy = malloc(nbr_atoms * sizeof (double));
    double *fz = malloc(nbr_atoms * sizeof (double));

    double *density = malloc(nbr_atoms * sizeof (double));
    double *dUembed_drho = malloc(nbr_atoms * sizeof (double));

    rcut = 6.06;
    rcut_sq = rcut * rcut;

    cell_length_inv = 1 / cell_length;
    cell_length_sq = cell_length * cell_length;

    for (i = 0; i < nbr_atoms; i++){
    sx[i] = positions[i][0] * cell_length_inv;
    sy[i] = positions[i][1] * cell_length_inv;
    sz[i] = positions[i][2] * cell_length_inv;
    }

    //每个原子初始电子密度和力
    for (i = 0; i < nbr_atoms; i++){
    density[i] = 0;
    fx[i] = 0;
    fy[i] = 0;
    fz[i] = 0;
    }


    //第i个原子周期性转换
    for (i = 0; i < nbr_atoms; i++) {
    /* Periodically translate coords of current particle to positive quadrants */
    sxi = sx[i] - floor(sx[i]); //floor 向下取整,2.8 -> 2.0
    syi = sy[i] - floor(sy[i]);
    szi = sz[i] - floor(sz[i]);

    densityi = density[i];

    /* Loop over other atoms. */
    for (j = i + 1; j < nbr_atoms; j++) {
    /* Periodically translate atom j to positive quadrants and calculate distance to it. */
    //第j个原子周期性处理,以及计算原子i和j的距离
    sxij = sxi - (sx[j] - floor(sx[j]));
    syij = syi - (sy[j] - floor(sy[j]));
    szij = szi - (sz[j] - floor(sz[j]));

    /* Periodic boundary conditions. */
    //原子距离周期性处理
    sxij = sxij - (int)floor(sxij + 0.5);
    syij = syij - (int)floor(syij + 0.5);
    szij = szij - (int)floor(szij + 0.5);

    /* squared distance between atom i and j */
    //i,j原子距离的平方和
    rij_sq = cell_length_sq * (sxij*sxij + syij*syij + szij*szij);

    /* Add force and energy contribution if distance between atoms smaller than rcut */
    //i,j原子的距离小于截断平方,则添加力和能量贡献
    if (rij_sq < rcut_sq) {
    rij = sqrt(rij_sq);
    //样条评估,返回密度
    dens = splineEval(rij, electron_density, ELECTRON_DENSITY_ROWS);
    //原子i和j的密度
    densityi += dens;
    density[j] += dens;
    }
    }
    density[i] = densityi;
    }

    /* Loop over atoms to calculate derivative of embedding function
    and embedding function. */
    //计算每个原子的嵌入函数导数和嵌入函数
    for (i = 0; i < nbr_atoms; i++) {
    dUembed_drho[i] = splineEvalDiff(density[i], embedding_energy, EMBEDDING_ENERGY_ROWS);
    }

    /* Compute forces on atoms. */
    /* Loop over atoms again :-(. */

    //计算每个原子的力
    for (i = 0; i < nbr_atoms; i++) {
    /* Periodically translate coords of current particle to positive quadrants */
    //周期性转换,同上
    sxi = sx[i] - floor(sx[i]);
    syi = sy[i] - floor(sy[i]);
    szi = sz[i] - floor(sz[i]);

    densityi = density[i];

    /* Loop over other atoms. */
    for (j = i + 1; j < nbr_atoms; j++) {
    /* Periodically translate atom j to positive quadrants and calculate distance to it. */
    //同上
    sxij = sxi - (sx[j] - floor(sx[j]));
    syij = syi - (sy[j] - floor(sy[j]));
    szij = szi - (sz[j] - floor(sz[j]));

    /* Periodic boundary conditions. */
    sxij = sxij - (int)floor(sxij + 0.5);
    syij = syij - (int)floor(syij + 0.5);
    szij = szij - (int)floor(szij + 0.5);

    /* squared distance between atom i and j */
    rij_sq = cell_length_sq * (sxij*sxij + syij*syij + szij*szij);

    /* Add force and energy contribution if distance between atoms smaller than rcut */
    //同上
    if (rij_sq < rcut_sq) {
    rij = sqrt(rij_sq);
    dUpair_dr = splineEvalDiff(rij, pair_potential, PAIR_POTENTIAL_ROWS);
    drho_dr = splineEvalDiff(rij, electron_density, ELECTRON_DENSITY_ROWS);

    /* Add force contribution from i-j interaction */
    //计算i,j原子之间的力,单位平均力=能量/i,j原子距离
    force = -(dUpair_dr + (dUembed_drho[i] + dUembed_drho[j])*drho_dr) / rij;
    //各个坐标方向,i和j原子的所受的力,这里是方形晶胞,因此各个方向力的大小相同
    fx[i] += force * sxij * cell_length;
    fy[i] += force * syij * cell_length;
    fz[i] += force * szij * cell_length;

    fx[j] -= force * sxij * cell_length;
    fy[j] -= force * syij * cell_length;
    fz[j] -= force * szij * cell_length;
    }
    }
    }

    //返回第i个原子的受力,x,y,z
    for (i = 0; i < nbr_atoms; i++){
    forces[i][0] = fx[i];
    forces[i][1] = fy[i];
    forces[i][2] = fz[i];
    }

    //释放内存
    free(sx); free(sy); free(sz); sx = NULL; sy = NULL; sz = NULL;
    free(fx); free(fy); free(fz); fx = NULL; fy = NULL; fz = NULL;
    free(density); density = NULL;
    free(dUembed_drho); dUembed_drho = NULL;

    }

    /* Returns the potential energy */
    //得到势能
    double get_energy_AL(double positions[][3], double cell_length, int nbr_atoms)
    {
    int i, j;
    double cell_length_inv, cell_length_sq;
    double rcut, rcut_sq;
    double energy;
    double densityi, dens;
    double sxi, syi, szi, sxij, syij, szij, rij, rij_sq;

    double *sx = malloc(nbr_atoms * sizeof (double));
    double *sy = malloc(nbr_atoms * sizeof (double));
    double *sz = malloc(nbr_atoms * sizeof (double));

    double *density = malloc(nbr_atoms * sizeof (double));

    rcut = 6.06;
    rcut_sq = rcut * rcut;

    cell_length_inv = 1 / cell_length;
    cell_length_sq = cell_length * cell_length;

    for (i = 0; i < nbr_atoms; i++){
    sx[i] = positions[i][0] * cell_length_inv;
    sy[i] = positions[i][1] * cell_length_inv;
    sz[i] = positions[i][2] * cell_length_inv;
    }

    for (i = 0; i < nbr_atoms; i++){
    density[i] = 0;
    }

    //初始能量0
    energy = 0;

    for (i = 0; i < nbr_atoms; i++) {
    /* Periodically translate coords of current particle to positive quadrants */
    sxi = sx[i] - floor(sx[i]);
    syi = sy[i] - floor(sy[i]);
    szi = sz[i] - floor(sz[i]);

    densityi = density[i];

    /* Loop over other atoms. */
    for (j = i + 1; j < nbr_atoms; j++) {
    /* Periodically translate atom j to positive quadrants and calculate distance to it. */
    sxij = sxi - (sx[j] - floor(sx[j]));
    syij = syi - (sy[j] - floor(sy[j]));
    szij = szi - (sz[j] - floor(sz[j]));

    /* Periodic boundary conditions. */
    sxij = sxij - (int)floor(sxij + 0.5);
    syij = syij - (int)floor(syij + 0.5);
    szij = szij - (int)floor(szij + 0.5);

    /* squared distance between atom i and j */
    rij_sq = cell_length_sq * (sxij*sxij + syij*syij + szij*szij);

    /* Add force and energy contribution if distance between atoms smaller than rcut */
    if (rij_sq < rcut_sq) {
    rij = sqrt(rij_sq);
    dens = splineEval(rij, electron_density, ELECTRON_DENSITY_ROWS);
    densityi += dens;
    density[j] += dens;

    /* Add energy contribution from i-j interaction */
    energy += splineEval(rij, pair_potential, PAIR_POTENTIAL_ROWS);

    }
    }
    density[i] = densityi;
    }

    /* Loop over atoms to calculate derivative of embedding function
    and embedding function. */
    //能量求和
    for (i = 0; i < nbr_atoms; i++) {
    energy += splineEval(density[i], embedding_energy, EMBEDDING_ENERGY_ROWS);
    }

    //释放内存
    free(sx); free(sy); free(sz); sx = NULL; sy = NULL; sz = NULL;
    free(density); density = NULL;

    //返回总能量
    return(energy);

    }

    /* Returns the virial */
    //得到维里
    double get_virial_AL(double positions[][3], double cell_length, int nbr_atoms)
    {
    int i, j;
    double cell_length_inv, cell_length_sq;
    double rcut, rcut_sq;
    double virial;
    double densityi, dens, drho_dr, force;
    double dUpair_dr;
    double sxi, syi, szi, sxij, syij, szij, rij, rij_sq;

    double *sx = malloc(nbr_atoms * sizeof (double));
    double *sy = malloc(nbr_atoms * sizeof (double));
    double *sz = malloc(nbr_atoms * sizeof (double));

    double *density = malloc(nbr_atoms * sizeof (double));
    double *dUembed_drho = malloc(nbr_atoms * sizeof (double));

    rcut = 6.06;
    rcut_sq = rcut * rcut;

    cell_length_inv = 1 / cell_length;
    cell_length_sq = cell_length * cell_length;

    for (i = 0; i < nbr_atoms; i++){
    sx[i] = positions[i][0] * cell_length_inv;
    sy[i] = positions[i][1] * cell_length_inv;
    sz[i] = positions[i][2] * cell_length_inv;
    }

    for (i = 0; i < nbr_atoms; i++){
    density[i] = 0;
    }

    for (i = 0; i < nbr_atoms; i++) {
    /* Periodically translate coords of current particle to positive quadrants */
    sxi = sx[i] - floor(sx[i]);
    syi = sy[i] - floor(sy[i]);
    szi = sz[i] - floor(sz[i]);

    densityi = density[i];

    /* Loop over other atoms. */
    for (j = i + 1; j < nbr_atoms; j++) {
    /* Periodically translate atom j to positive quadrants and calculate distance to it. */
    sxij = sxi - (sx[j] - floor(sx[j]));
    syij = syi - (sy[j] - floor(sy[j]));
    szij = szi - (sz[j] - floor(sz[j]));

    /* Periodic boundary conditions. */
    sxij = sxij - (int)floor(sxij + 0.5);
    syij = syij - (int)floor(syij + 0.5);
    szij = szij - (int)floor(szij + 0.5);

    /* squared distance between atom i and j */
    rij_sq = cell_length_sq * (sxij*sxij + syij*syij + szij*szij);

    /* Add force and energy contribution if distance between atoms smaller than rcut */
    if (rij_sq < rcut_sq) {
    rij = sqrt(rij_sq);
    dens = splineEval(rij, electron_density, ELECTRON_DENSITY_ROWS);
    densityi += dens;
    density[j] += dens;
    }
    }
    density[i] = densityi;
    }

    /* Loop over atoms to calculate derivative of embedding function
    and embedding function. */
    for (i = 0; i < nbr_atoms; i++) {
    dUembed_drho[i] = splineEvalDiff(density[i], embedding_energy, EMBEDDING_ENERGY_ROWS);
    }

    /* Compute forces on atoms. */
    /* Loop over atoms again :-(. */

    //初始维里为0
    virial = 0;

    for (i = 0; i < nbr_atoms; i++) {
    /* Periodically translate coords of current particle to positive quadrants */
    sxi = sx[i] - floor(sx[i]);
    syi = sy[i] - floor(sy[i]);
    szi = sz[i] - floor(sz[i]);

    densityi = density[i];

    /* Loop over other atoms. */
    for (j = i + 1; j < nbr_atoms; j++) {
    /* Periodically translate atom j to positive quadrants and calculate distance to it. */
    sxij = sxi - (sx[j] - floor(sx[j]));
    syij = syi - (sy[j] - floor(sy[j]));
    szij = szi - (sz[j] - floor(sz[j]));

    /* Periodic boundary conditions. */
    sxij = sxij - (int)floor(sxij + 0.5);
    syij = syij - (int)floor(syij + 0.5);
    szij = szij - (int)floor(szij + 0.5);

    /* squared distance between atom i and j */
    rij_sq = cell_length_sq * (sxij*sxij + syij*syij + szij*szij);

    /* Add force and energy contribution if distance between atoms smaller than rcut */
    if (rij_sq < rcut_sq) {
    rij = sqrt(rij_sq);
    dUpair_dr = splineEvalDiff(rij, pair_potential, PAIR_POTENTIAL_ROWS);
    drho_dr = splineEvalDiff(rij, electron_density, ELECTRON_DENSITY_ROWS);

    /* Add virial contribution from i-j interaction */
    force = -(dUpair_dr + (dUembed_drho[i] + dUembed_drho[j])*drho_dr) / rij;

    //维里累加,维里=力x距离平方
    virial += force * rij_sq;
    }
    }
    }

    //维里/3
    virial /= 3.0;

    //释放内存
    free(sx); free(sy); free(sz); sx = NULL; sy = NULL; sz = NULL;
    free(density); density = NULL;
    free(dUembed_drho); dUembed_drho = NULL;

    return(virial);

    }
  • MD_main.c文件,主函数,得到体系势能与体积的变化关系。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    /*
    MD_main.c

    Created by Anders Lindman on 2013-10-31.
    */

    #include <stdio.h>
    #include <math.h>
    #include <stdlib.h>
    #include <time.h>
    #include "initfcc.h"
    #include "alpotential.h"

    /* Main program */
    int main()
    {

    /*
    Code for generating a uniform random number between 0 and 1. srand should only
    be called once.
    */
    /*
    srand(time(NULL));
    double random_value;
    random_value = (double) rand() / (double) RAND_MAX;
    */

    /*
    Descriptions of the different functions in the files initfcc.c and
    alpotential.c are listed below.
    */

    /*
    Function that generates a fcc lattice in units of [Å]. Nc is the number of
    primitive cells in each direction and a0 is the lattice parameter. The
    positions of all the atoms are stored in pos which should be a matrix of the
    size N x 3, where N is the number of atoms. The first, second and third column
    correspond to the x,y and z coordinate respectively.
    */
    FILE *file;
    file = fopen("plot.dat", "w");

    for (int i = 0; i < 12; i++)
    {

    int Nc = 4;
    int N = 4 * Nc * Nc * Nc;
    double v0 = 63 + i / 2.0;
    double a0 = pow(v0, 1.0 / 3.0);
    double pos[N][3];

    init_fcc(pos, Nc, a0);

    /*
    Function that calculates the potential energy in units of [eV]. pos should be
    a matrix containing the positions of all the atoms, L is the length of the
    supercell and N is the number of atoms.
    */

    double energy;
    energy = get_energy_AL(pos, Nc * a0, N);
    fprintf(file, "%.4f %.4f\n", a0 * a0 * a0, energy / N * 4);
    }
    fclose(file);
    /*
    Function that calculates the virial in units of [eV]. pos should be a matrix
    containing the positions of all the atoms, L is the length of the supercell
    and N is the number of atoms.
    */
    /*
    double virial;
    virial = get_virial_AL(pos, L, N);
    */

    /*
    Function that calculates the forces on all atoms in units of [eV/Å]. the
    forces are stored in f which should be a matrix of size N x 3, where N is the
    number of atoms and column 1,2 and 3 correspond to the x,y and z component of
    the force resepctively . pos should be a matrix containing the positions of
    all the atoms, L is the length of the supercell and N is the number of atoms.
    */
    /*
    get_forces_AL(f,pos, L, N);
    */
    }

编译

Linux下:

1
gcc -O3 -Wall initfcc.c alpotential.c MD_main.c -o md -lm

或者Windows下

1
icl /O3 initfcc.c alpotential.c MD_main.c -o md

得到md或者md.exe可执行文件。

运行及数据绘图

Linux下直接:

1
./md

Windows下,打开CMD,执行:

1
md.exe

或者直接双击也行。
生成数据文件plot.dat


绘图得到势能与体积的关系如下所示:

讨论

如上图所示,当单位晶胞体积为$V_0=65.6 \dot{A}^3$时势能最小,此时原子的吸引力和排斥力相互抵消。这对应于$T=0K时的$晶格常数a0,也就是$(V_0)^\frac{1}{3}=4.03\dot{A}$。测量表明a0的真实值为$4.0317\dot{A}$。