クーロンの法則力算出使ってCコード

Question

私は10個のビーズを含んbeads.Both鎖の2つの鎖の間の対相互作用（クーロン力）を計算するためにC言語で簡単なプログラムを記述しようとしていますが、1つの鎖は、正に帯電し、他方はです負に帯電した。私はx、y、z座標と2つのファイルから電荷情報を読み込んでいます（atom.psfの7番目の列からの電荷とbeads.pdbの6番目の7番目と8番目の列からの座標）。ここ は、ファイルの内容は以下のとおりです。

ATOM  1 A CHN  1  1.000 0.000 0.000 1.00 0.00  A 
ATOM  2 A CHN  1  2.000 0.000 0.000 1.00 0.00  A 
ATOM  3 A CHN  1  3.000 0.000 0.000 1.00 0.00  A 
ATOM  4 A CHN  1  4.000 0.000 0.000 1.00 0.00  A 
ATOM  5 A CHN  1  5.000 0.000 0.000 1.00 0.00  A 
ATOM  6 A CHN  1  6.000 0.000 0.000 1.00 0.00  A 
ATOM  7 A CHN  1  7.000 0.000 0.000 1.00 0.00  A 
ATOM  8 A CHN  1  8.000 0.000 0.000 1.00 0.00  A 
ATOM  9 A CHN  1  9.000 0.000 0.000 1.00 0.00  A 
ATOM  10 A CHN  1  10.000 0.000 0.000 1.00 0.00  A 
ATOM  11 A CHN  2  1.000 80.000 0.000 1.00 0.00  A 
ATOM  12 A CHN  2  2.000 80.000 0.000 1.00 0.00  A 
ATOM  13 A CHN  2  3.000 80.000 0.000 1.00 0.00  A 
ATOM  14 A CHN  2  4.000 80.000 0.000 1.00 0.00  A 
ATOM  15 A CHN  2  5.000 80.000 0.000 1.00 0.00  A 
ATOM  16 A CHN  2  6.000 80.000 0.000 1.00 0.00  A 
ATOM  17 A CHN  2  7.000 80.000 0.000 1.00 0.00  A 
ATOM  18 A CHN  2  8.000 80.000 0.000 1.00 0.00  A 
ATOM  19 A CHN  2  9.000 80.000 0.000 1.00 0.00  A 
ATOM  20 A CHN  2  10.000 80.000 0.000 1.00 0.00  A 

beads.pdb

20 !NATOM 
     1 A 1 CHN A  A 1.000000  0.0000   0 
     2 A 1 CHN A  A 1.000000  0.0000   0 
     3 A 1 CHN A  A 1.000000  0.0000   0 
     4 A 1 CHN A  A 1.000000  0.0000   0 
     5 A 1 CHN A  A 1.000000  0.0000   0 
     6 A 1 CHN A  A 1.000000  0.0000   0 
     7 A 1 CHN A  A 1.000000  0.0000   0 
     8 A 1 CHN A  A 1.000000  0.0000   0 
     9 A 1 CHN A  A 1.000000  0.0000   0 
     10 A 1 CHN A  A 1.000000  0.0000   0 
     11 A 2 CHN A  A -1.000000  0.0000   0 
     12 A 2 CHN A  A -1.000000  0.0000   0 
     13 A 2 CHN A  A -1.000000  0.0000   0 
     14 A 2 CHN A  A -1.000000  0.0000   0 
     15 A 2 CHN A  A -1.000000  0.0000   0 
     16 A 2 CHN A  A -1.000000  0.0000   0 
     17 A 2 CHN A  A -1.000000  0.0000   0 
     18 A 2 CHN A  A -1.000000  0.0000   0 
     19 A 2 CHN A  A -1.000000  0.0000   0 
     20 A 2 CHN A  A -1.000000  0.0000   0 

     18 !NBOND: bonds 
     1  2  2  3  3  4  4  5 
     5  6  6  7  7  8  8  9 
     9  10  11  12  12  13  13  14 
     14  15  15  16  16  17  17  18 
     18  19  19  20 

atom.psfは、私は私の最終的な出力とのトラブルを抱えています。すべてのタイムステップ（t = 0〜100）で、私は20の原子の座標を書く必要があります。私の試用コードは以下の通りです。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <string.h> 

#define Epsilon0 8.85e-12 // Permittivity of free space (C^2/(N m^2)) 
#define Constant (1/(4*M_PI*Epsilon0)) // Useful constant 
#define gamma 100.0 
#define ROW 20 

int row, col, j, t; 
float x[20], y[20], z[20], q[20], dx, dy, dz, R, Fx, Fy, Fz, F, v, shift; 

int main() 
{ 
    FILE *psf=fopen("atom.psf", "r"); 
    FILE *pdb=fopen("beads.pdb", "r"); 
    FILE *fout=fopen("out.txt", "a"); 
    FILE *fout2=fopen("coord.dump", "a"); 
    fprintf(fout2, "ITEM: TIMESTEP\n 100\n"); 
    fprintf(fout2, "ITEM: NUMBER OF ATOMS\n 20\n"); 
    fprintf(fout2, "ITEM: ATOMS id type x y z\n"); 
    char buffer[1024]; 

    fgets(buffer, 1024, psf); 

    int i = 0; 

    for(i=0 ; (i<ROW) && (psf != NULL); i++) 
    { 
    fscanf (psf,"%*8d%*4s%*4d%*6s%*5s%*6s%11f%*14f%*12d", &q[i]); 
    fscanf (pdb,"%*4s%*7d%*5s%*4s%*6d%12f%8f%8f%*6f%*6f%*9s", &x[i], &y[i], &z[i]); 
    } 

    for (t=0; t<100; t++) 
    { 
     //F = 0.0; 
     v = F/gamma; 
     shift = v*t; 
     x[i] = x[i] + shift; 
     y[i] = y[i] + shift; 
     z[i] = z[i] + shift; 
     for(i=0; i<ROW; i++) 
     { 
     Fx = Fy = Fz = F = 0.0; 
     // Loop over other charges to compute force on this charge 
     for (j=0 ; j<ROW ; j++) 
     { 
     //simply skip this itearation 
      if(i == j) 
      continue; 
     // Compute the components of vector distance between two charges 
      dx = x[i] - x[j]; 
      dy = y[i] - y[j]; 
      dz = z[i] - z[j]; 
      R = sqrt(dx*dx + dy*dy + dz*dz); 

      // Compute the x and y components of the force between 
      // these two charges using Coulomb's law 
      Fx += Constant*q[i]*q[j]*dx/(R*R*R); 
      Fy += Constant*q[i]*q[j]*dy/(R*R*R); 
      Fz += Constant*q[i]*q[j]*dz/(R*R*R); 
     } 

      F = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); 

    fprintf(fout, "%d %3.3g %3.3g %3.3g %3.3g\n", i+1, Fx, Fy, Fz, F); 
    //fprintf(fout2, "%d %3.3g %3.3g %3.3g\n", i+1, x[i], y[i], z[i]); 
    } 
    fprintf(fout2, "%d %3.3g %3.3g %3.3g\n", i, x[i], y[i], z[i]); 
    } 
} 

Answer 1

を完了したときF[i]を計算し、その後、すなわちFx[i] +=、各interationためaccumlateする変更を検討：

• イプシロンと定数として宣言する必要があります定数： "の#define Epsilon0の8.85e-12" "の#define定数（1 /（4 *をPI * Epsilon0））"」

• ループ内で "i！= j"を移動します。その状態になるとループが止まるのを防ぎます。完全に停止するのではなく、その反復をスキップするだけです。

• は、あなたは本当にXI、YI、zij、Fxの、FY、Fzを、Rijのアレイを必要としません。代わりに単純変数をプレースホルダとして使用できます。内側のループの各反復の間に

• は、あなたは、ビーズjへの部分的な力デュを計算することになります。その力を累積力に加える必要があります。 あなたは、ループの外でのFx、FyとFzとを宣言することができることをやるためには、最初のループ内で0に3変数を初期化して、内側のループ内で部分的に力を加えます。

• 移動内側ループの外側関数fprintf（ただし、外部ループ内に保つ）

更新が

このコードは、Iを処理しない/ Oエラー

#include <stdio.h> 
    #include <stdlib.h> 
    #include <math.h> 
    #include <string.h> 

    #define Epsilon0 8.85e-12 // Permittivity of free space (C^2/(N m^2)) 
    #define Constant (1/(4*M_PI*Epsilon0)) // Useful constant 

    #define ROW 20 

    int row, col, j; 
    float x[20], y[20], z[20], q[20], dx, dy, dz, R, Fx, Fy, Fz, F; 

    int main() 
    { 
     FILE *psf=fopen("atom.psf", "r"); 
     FILE *pdb=fopen("beads.pdb", "r"); 
     FILE *fout=fopen("out.txt", "w"); 
     char buffer[1024]; 

     fgets(buffer, 1024, psf); 

     int i = 0; 

     for(; (i<ROW) && (psf != NULL); i++) 
     { 
     fscanf (psf,"%*8d%*4s%*4d%*6s%*5s%*6s%11f%*14f%*12d", &q[i]); 
     fscanf (pdb,"%*4s%*7d%*5s%*4s%*6d%12f%8f%8f%*6f%*6f%*9s", &x[i], &y[i], &z[i]); 
     } 

     for(i=0; i<ROW; i++) 
     { 
     Fx = Fy = Fz = F = 0.0; 
     // Loop over other charges to compute force on this charge 
     for (j=0 ; j<ROW ; j++) 
     { 
      //simply skip this itearation 
      if(i == j) 
      continue; 
      // Compute the components of vector distance between two charges 
      dx = x[i] - x[j]; 
      dy = y[i] - y[j]; 
      dz = z[i] - z[j]; 
      R = sqrt(dx*dx + dy*dy + dz*dz); 

      // Compute the x and y components of the force between 
      // these two charges using Coulomb's law 
      Fx += Constant*q[i]*q[j]*dx/(R*R*R); 
      Fy += Constant*q[i]*q[j]*dy/(R*R*R); 
      Fz += Constant*q[i]*q[j]*dz/(R*R*R); 
     } 

     F = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); 

     fprintf(fout, "%d %g %g %g, %g\n", i+1, Fx, Fy, Fz, F); 
     } 
    } 

Answer 2

ループ条件としてi == jがfalseの場合

`for (j = 0; j<ROW && i != j; j++)` 

が終了します、forループ

ザ・を助けるかもしれない何かを。私はあるときに充電電流をスキップするループのための、すなわち内側のif文i == j

for (j = 0; j<ROW ; j++) 
    { 
     if (i != j) 
     { 
      DO CALCULATION... ` 

または、と思うときjは0

使用する場合例えば、1は、一回のループを繰り返し処理しますFx[i]は今まで、それは各反復の割り当てではなく和であるように、単一のj電荷に起因する力を格納しています。

ここで私はあなたのコードに変更します何このforループが

Answer 3

これはおそらく、さまざまな情報を含む多数の別々の配列を追跡しようとする代わりに、li feは、それぞれビーズの必要な情報をキャプチャし、という単一の配列を作成すると、structを作成するだけで簡単に作成できます。たとえば、各ビードに作用するx、y、zの位置、電荷、および成分の力をキャプチャしたいとします。 （私は合計力も含めたが、これはオプションである）。あなたは、単にbeadsの配列（あなたが情報を必要とする各ビーズに1つ）を作成し、あなたのコードで

typedef struct { 
    float x, y, z, chrg, fx, fy, fz, f; 
} beads; 

：各ビーズのためのあなたのstructは（私はbeadsそれを呼ばれる）のような単純なものでした。例えば、20要素を持つビーズの配列を作成するには、次のような何かができる：

#define Eo 8.85e-12F 
#define KEair (1/(4*M_PI*Eo)) 

enum { NATM = 20, MAXL = 128 }; /* constants for number of beads/line len */ 
... 
    beads atoms[NATM] = {{ .x = 0.0 }}; 

は、今、私たちは、構造体の配列は、情報を格納するatomsを求めているあなたは、あなたのファイルのそれぞれから情報を読み取ることができます。各ビードのx, y, zの位置とchargeの位置を保存します。その後、他のすべてのビードによる力を計算し、残りの力成分と各構造体の合計メンバーの情報を合計することによって、各ビードに作用する力を計算することができます。

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 

#define Eo 8.85e-12F 
#define KEair (1/(4*M_PI*Eo)) 

enum { NATM = 20, MAXL = 128 }; 

typedef struct { 
    float x, y, z, chrg, fx, fy, fz, f; 
} beads; 

FILE *xfopen (const char *fn, const char *mode); 
float sqrt_fisr (float x); 

int main (int argc, char **argv) { 

    beads atoms[NATM] = {{ .x = 0.0 }}; 
    size_t i; 
    char buf[MAXL] = ""; /* open/read atom.psf */ 
    FILE *fp = xfopen (argc > 1 ? argv[1] : "dat/atom.psf", "r"); 

    fgets (buf, MAXL, fp); /* read/discard 1st line */ 
    for (i = 0; i < NATM && fgets (buf, MAXL, fp); i++) { /* read/parse data */ 
     if (sscanf (buf, "%*s %*s %*s %*s %*s %*s %f", &atoms[i].chrg) != 1) { 
      fprintf (stderr, "error: read of charge failed, atom[%zu].\n", i); 
      return 1; 
     } 
    } 
    fclose (fp); 
    if (i != NATM) { /* validate NATM lines read */ 
     fprintf (stderr, "error: only '%zu' charge values read.\n", i); 
     return 1; 
    } 
     /* open/read beads.pdb */ 
    fp = xfopen (argc > 2 ? argv[2] : "dat/beads.pdb", "r"); 

    for (i = 0; i < NATM && fgets (buf, MAXL, fp); i++) { /* read/parse data */ 
     if (sscanf (buf, "%*s %*s %*s %*s %*s %f %f %f", &atoms[i].x, 
      &atoms[i].y, &atoms[i].z) != 3) { 
      fprintf (stderr, "error: read of position failed, atom[%zu].\n", i); 
      return 1; 
     } 
    } 
    fclose (fp); 
    if (i != NATM) { /* validate NATM lines read */ 
     fprintf (stderr, "error: only '%zu' position values read.\n", i); 
     return 1; 
    } 

    for (i = 0; i < NATM; i++) {    /* for each bead */ 
     for (size_t j = 0; j < NATM; j++) {  /* compute force from every other */ 
      if (i == j) continue;    /* excluding itself */ 
      float dx = atoms[j].x - atoms[i].x, /* calculate component distances */ 
        dy = atoms[j].y - atoms[i].y, 
        dz = atoms[j].z - atoms[i].z, 
        d = sqrt (dx * dx + dy * dy + dz * dz); /* total distance */ 

      /* compute component and total forces acting on each bead (sum) */ 
      atoms[i].fx += (KEair * atoms[i].chrg *atoms[j].chrg * dx)/(d * d * d); 
      atoms[i].fy += (KEair * atoms[i].chrg *atoms[j].chrg * dy)/(d * d * d); 
      atoms[i].fz += (KEair * atoms[i].chrg *atoms[j].chrg * dz)/(d * d * d); 
      atoms[i].f += (KEair * atoms[i].chrg *atoms[j].chrg)/(d * d); 
     } 
    } 

    for (i = 0; i < NATM; i++) /* output forces on each bead (component and total) */ 
     printf (" atom[%2zu] %5.2f %5.2f %5.2f %+.2f %15.2f %15.2f %5.2f %15.2f\n", 
       i, atoms[i].x, atoms[i].y, atoms[i].z, atoms[i].chrg, 
       atoms[i].fx, atoms[i].fy, atoms[i].fz, atoms[i].f); 

    return 0; 
} 

/** simple fopen with error check */ 
FILE *xfopen (const char *fn, const char *mode) 
{ 
    FILE *fp = fopen (fn, mode); 

    if (!fp) { 
     fprintf (stderr, "xfopen() error: file open failed '%s'.\n", fn); 
     exit (EXIT_FAILURE); 
    } 

    return fp; 
} 

例：あなたは、次のような何かを行うことができ、

for (i = 0; i < NATM; i++) {    /* for each bead */ 
     for (size_t j = 0; j < NATM; j++) {  /* compute force from every other */ 
      if (i == j) continue;    /* excluding itself */ 
      float dx = atoms[j].x - atoms[i].x, /* calculate component distances */ 
        dy = atoms[j].y - atoms[i].y, 
        dz = atoms[j].z - atoms[i].z, 
        d = sqrt (dx * dx + dy * dy + dz * dz); /* total distance */ 

      /* compute component and total forces acting on each bead (sum) */ 
      atoms[i].fx += (KEair * atoms[i].chrg *atoms[j].chrg * dx)/(d * d * d); 
      atoms[i].fy += (KEair * atoms[i].chrg *atoms[j].chrg * dy)/(d * d * d); 
      atoms[i].fz += (KEair * atoms[i].chrg *atoms[j].chrg * dz)/(d * d * d); 
      atoms[i].f += (KEair * atoms[i].chrg *atoms[j].chrg)/(d * d); 
     } 
    } 

（あなたがSystem of discrete chargesで合計へのアプローチを確認することができます）

完全にそれを置く：あなたは何かのように行うことができます使用/出力

$ ./bin/coulomb 
atom[ 0] 1.00 0.00 0.00 +1.00 13844509696.00  -13956823.00 0.00 13831302144.00 
atom[ 1] 2.00 0.00 0.00 +1.00  4741866496.00  -13982750.00 0.00 22712082432.00 
atom[ 2] 3.00 0.00 0.00 +1.00  2353591552.00  -14002249.00 0.00 24819521536.00 
atom[ 3] 4.00 0.00 0.00 +1.00  1171170304.00  -14015272.00 0.00 25635096576.00 
atom[ 4] 5.00 0.00 0.00 +1.00  359584800.00  -14021791.00 0.00 25947308032.00 
atom[ 5] 6.00 0.00 0.00 +1.00  -359584608.00  -14021791.00 0.00 25947308032.00 
atom[ 6] 7.00 0.00 0.00 +1.00 -1171170944.00  -14015272.00 0.00 25635096576.00 
atom[ 7] 8.00 0.00 0.00 +1.00 -2353592320.00  -14002249.00 0.00 24819521536.00 
atom[ 8] 9.00 0.00 0.00 +1.00 -4741866496.00  -13982751.00 0.00 22712082432.00 
atom[ 9] 10.00 0.00 0.00 +1.00 -13844510720.00  -13956824.00 0.00 13831303168.00 
atom[10] 1.00 80.00 0.00 -1.00 13844507648.00  13956823.00 0.00 13831302144.00 
atom[11] 2.00 80.00 0.00 -1.00  4741867008.00  13982750.00 0.00 22712080384.00 
atom[12] 3.00 80.00 0.00 -1.00  2353592064.00  14002249.00 0.00 24819521536.00 
atom[13] 4.00 80.00 0.00 -1.00  1171170944.00  14015272.00 0.00 25635096576.00 
atom[14] 5.00 80.00 0.00 -1.00  359585056.00  14021791.00 0.00 25947308032.00 
atom[15] 6.00 80.00 0.00 -1.00  -359584864.00  14021791.00 0.00 25947308032.00 
atom[16] 7.00 80.00 0.00 -1.00 -1171170560.00  14015272.00 0.00 25635096576.00 
atom[17] 8.00 80.00 0.00 -1.00 -2353591808.00  14002249.00 0.00 24819521536.00 
atom[18] 9.00 80.00 0.00 -1.00 -4741867008.00  13982751.00 0.00 22712080384.00 
atom[19] 10.00 80.00 0.00 -1.00 -13844509696.00  13956824.00 0.00 13831304192.00 

出力を視覚的に見ると、力は各ストリングまたはビーズの中心の端から最大まで増加します。これはチェックします。ビーズに作用するX方向の力は両端で最大であり、中点で - 小切手でもあります。たとえば、次のように

全磁力演技X方向に作用する力

各ビーズ上のルック物事を超えるとさせて頂いて

​​

各ビーズにあなたは何か質問がある場合は知っている。