nabaneeta mukhopadhyay - 10 months ago 117
C Question

# coulomb's law force calculation using c code

I am trying to write a simple program in C to calculate pairwise interaction (Coulomb force) between two chain of beads.Both chains contain 10 beads but one chain is positively charged and the other one is negatively charged. I am reading the x,y,z coordinates and the charge information from two files (Charge from 7th column of atom.psf and coordinates from 6th 7th and 8th columns of beads.pdb).
Here are the contents of the files:
atom.psf

``````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      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
``````

I am having trouble with my final output. There should be a list of 20 force values (each bead is experiencing force due to other 19 beads). I believe I am calculating individual forces but not sure how to sum them up for each bead. My trial code is given below.

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

const int ROW=20;

int row, col, j;
float x[20], y[20], z[20], q[20], xij[20], yij[20], zij[20], Rij[20], Fx[20], Fy[20], Fz[20], F[20];

float pi = 3.14;
float Epsilon0 = 8.85e-12;    // Permittivity of free space (C^2/(N m^2))

int main()
{
FILE *psf=fopen("atom.psf", "r");
FILE *fout=fopen("out.txt", "w");
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(i=0; i<ROW; i++)
{

// Loop over other charges to compute force on this charge

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

// Compute the components of vector distance between two charges
xij[i] = x[i] - x[j];
yij[i] = y[i] - y[j];
zij[i] = z[i] - z[j];
Rij[i] = sqrt(xij[i]*xij[i] + yij[i]*yij[i] +zij[i]*zij[i]);

float Constant = (1/(4*pi*Epsilon0));  // Useful constant

// Compute the x and y components of the force between
// these two charges using Coulomb's law
Fx[i] = Constant*q[i]*q[j]*xij[i]/(Rij[i]*Rij[i]*Rij[i]);
Fy[i] = Constant*q[i]*q[j]*yij[i]/(Rij[i]*Rij[i]*Rij[i]);
Fz[i] = Constant*q[i]*q[j]*zij[i]/(Rij[i]*Rij[i]*Rij[i]);

// Here I need to sum up the forces on each bead due to other 19 beads and I need help.
//  F[i] = Fx[i] + Fy[i] + Fz[i];

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

}
}
}
``````

This is probably a situation where instead of attempting to keep track of numerous separate arrays containing various information, life can be made a lot easier if you simply create a `struct` that captures the needed information for each bead and then create a single array of struct. For example, you want to capture the x,y,z position, charge, and component forces acting on each bead. (I included the total force as well, but that is optional). Your `struct` for each bead (I called it `beads`) could be as simple as:

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

In your code, you simply create an array of `beads` (one for each bead you need information on). For example, creating an array of beads with `20` elements, you could do something like the following:

``````#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 }};
``````

Now we have an array of struct called `atoms` to store information in. You can read information from each of your files and store the `x, y, z` positions for each bead as well as the `charge`. Then you can compute the force acting on each bead by computing the force due to every other bead and summing the information in the remaining force component and total members of each struct. You may do something like:

``````    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);
}
}
``````

(you can confirm the approach to the sum at System of discrete charges)

Putting it altogether, you could do something like the following:

``````#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;

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");

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;
}
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;
}
``````

Example Use/Output

``````\$ ./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
``````

Looking at the output graphically, the forces increase from the end to a maximum in the center of each string or beads -- that checks, and the forces acting upon the beads in the X-direction is the greatest at each end and reverses direction at midpoint -- also a check. For example:

Total Force Acting on Each Bead

Forces acting in the X direction on Each Bead

Look things over and let me know if you have any questions.

Source (Stackoverflow)