nabaneeta mukhopadhyay nabaneeta mukhopadhyay - 3 months ago 57
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


beads.pdb

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 *pdb=fopen("beads.pdb", "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]);



}
}
}

Answer

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

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

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

enter image description here

Forces acting in the X direction on Each Bead

enter image description here

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

Comments