Reputation: 47
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. At every timestep (t=0 to 100), I need to write the coordinates of 20 atoms. My trial code is given below.
#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]);
}
}
Upvotes: 0
Views: 3129
Answers (3)
Reputation: 26
Here's what I would change in your code:
Epsilon and Constant should be declared as constants: "#define Epsilon0 8.85e-12" "#define Constant (1/(4*pi*Epsilon0))""
Move the "i != j" inside the loop. That would prevent the loop from stopping if you hit that condition. You would just be skipping that iteration instead of stopping completely.
You don't really need the xi, yi , zij, Fx, Fy, Fz, Rij arrays. You can use simple variables as placeholders instead.
During each iteration of the inner loop, you would be calculating the partial force du to bead j. You need to add that force to the cumulative force. In order to do that you can declare Fx, Fy and Fz outside of the loops, initialize the 3 variable to 0 inside the first loop and then add partial forces inside the inner loop.
Move the fprintf outside of the inner loop (but keep it inside the external loop)
UPDATE
This code does not handle I/O errors
#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);
}
}
Upvotes: 1
Reputation: 84561
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
Forces acting in the X direction on Each Bead
Look things over and let me know if you have any questions.
Upvotes: 2
Reputation: 1
Something that may help
The for loop,
`for (j = 0; j<ROW && i != j; j++)`
Will exit when i == j as the loop condition is false. For example when i is 1 will only iterate through the loop once, when j is 0.
Use if statement inside the for loop to skip the current charge, i.e when i == j
for (j = 0; j<ROW ; j++)
{
if (i != j)
{
DO CALCULATION... `
Also, think that Fx[i] will only ever have store the force due to an single j charge as it is assignment for each iteration rather than a sum.
Consider changing to accumlate for each interation, i.e. Fx[i] += then calculate F[i] when this for loop has completed
Upvotes: 0
Related Questions
- Summing and storing variables of an array in a for loop C
- Coutning Table Using Summation and Inner for loop
- Sum of arrays in C
- Inputting value on 2D Float arrays with for loops
- Iterating through a multidimensional array and adding in floats
- How to Construct Multidimensional Array with Variable Assignment
- How to increment each element of 2-dimensional array by 2
- Performance issue in Finite difference method
- Using 2D Arrays in C
- Computationally efficient three dimensional arrays in C