个人中心
文章发布
退出
作者热门文章
- android - 多次调用 OnPrimaryClipChangedListener
- android - 无法更新 RecyclerView 中的 TextView 字段
- android.database.CursorIndexOutOfBoundsException : Index 0 requested, 光标大小为 0
- android - 使用 AppCompat 时，我们是否需要明确指定其 UI 组件(Spinner、EditText)颜色
24
4
我正在尝试用 C 编写一个简单的程序来计算两条珠链之间的成对相互作用(库仑力)。两条链都包含 10 个珠子，但一条链带正电，另一 strip 负电。我正在从两个文件中读取 x、y、z 坐标和电荷信息(电荷来自 atom.psf 的第 7 列，坐标来自 beads.pdb 的第 6、7 和 8 列)。以下是文件的内容:原子.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
珠子.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
我的最终输出有问题。在每个时间步(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]);
}
}
最佳答案
这可能是这样一种情况，与其尝试跟踪包含各种信息的大量单独数组，不如简单地创建一个 struct 来捕获每个数组所需的信息，这样生活就会轻松很多珠子，然后创建一个结构数组。例如，您想要捕获作用在每个珠子上的 x、y、z 位置、电荷和分力。 (我也包括总力，但这是可选的)。每个珠子的 struct(我称之为 beads)可以简单如下:
typedef struct {
float x, y, z, chrg, fx, fy, fz, f;
} beads;
在您的代码中，您只需创建一组 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。然后，您可以计算作用在每个珠子上的力，方法是计算每个其他珠子的力，并将剩余力分量和每个结构的总成员中的信息相加。你可以这样做:
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 确认求和的方法)
总而言之，您可以执行以下操作:
#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;
}
示例使用/输出
$ ./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 方向的力
检查一下，如果您有任何问题，请告诉我。
关于使用C代码计算库仑定律，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/38818455/
24
4
0
SQL 和一般开发的新手，我有一个表(COUNTRIES)，其中包含字段(INDEX、NAME、POPULATION、AREA) 通常我添加一个客户端(Delphi)计算字段(DENSITY)和 On
我想使用 calc(100%-100px)，但在我的 demo 中不起作用由于高度只接受像素，因此如何将此百分比值转换为像素。 最佳答案 以下将为您提供高度: $(window).height();
我正在尝试在 MySQL 中添加列并动态填充其他列。 例如我有一张表“数字”并具有第 1 列、第 2 列、第 3 列，这些总数应填充在第 4 列中 最佳答案 除非我误解了你的问题，否则你不只是在寻找:
我想返回简单计算的结果，但我不确定如何执行此操作。我的表格如下: SELECT COUNT(fb.engineer_id) AS `total_feedback`, SUM(fb.ra
我一直在尝试做这个程序，但我被卡住了，我仍然是一个初学者，任何帮助将不胜感激。我需要程序来做 打印一个 10 X 10 的表格，其中表格中的每个条目都是行号和列号的总和 包含一个累加器，用于计算所有表
这个计算背后一定有一些逻辑。但我无法得到它。普通数学不会导致这种行为。谁能帮我解释一下原因 printf ("float %f\n", 2/7 * 100.0); 结果打印 1.000000 为什么会
我想计算从 0 到 (n)^{1/2} - 1 的数字的 AND每个数字从 0 到 (n)^{1/2} - 1 .我想在 O(n) 中执行此操作时间，不能使用 XOR、OR、AND 运算。 具体来说，
如何在 Excel 中将公式放入自定义数字格式？例如(出于说明目的随机示例)， 假设我有以下数据: 输入 输出 在不编辑单元格中的实际数据的情况下，我想显示单元格中的值除以 2，并保留两位小数: 有没
每次我在 Flutter 应用程序中调用计算()时，我都会看到内存泄漏，据我所知，这基本上只是一种生成隔离的便捷方法。我的应用程序内存占用增加并且在 GC 之后永远不会减少。 我已将我的代码简化为仅调
我有数字特征观察 V1通过 V12用于目标变量 Wavelength .我想计算 Vx 之间的 RMSE列。数据格式如下。 每个变量“Vx”以 5 分钟的间隔进行测量。我想计算所有 Vx 变量的观测值
我正在寻找一种使用 C 语言计算文件中未知字符数的简单方法。谢谢你的帮助 最佳答案 POSIX 方式(可能是您想要的方式): off_t get_file_length( FILE *file ) {
我正在使用 Postgres，并且我正试图围绕如何在连续日期跨度中得出第一个开始日期的问题进行思考。例如 :- ID | Start Date | End Date =================
我有一个订单表格，我在其中使用 jQuery 计算插件来汇总总数。 此求和工作正常，但生成的“总和”存在问题。总之，我希望用逗号替换任何点。 代码的基础是； function ($this) {
我在使用 double 变量计算简单算术方程时遇到问题。 我有一个具有 double 属性 Value 的组件，我将此属性设置为 100。 然后我做一个简单的减法来检查这个值是否真的是 100: va
我在这里看到了一些关于 CRC 32 计算的其他问题。但没有一个让我满意，因此是这样。 openssl 库是否有任何用于计算 CRC32 的 api 支持？我已经在为 SHA1 使用 openssl，
当我在PHP日期计算中遇到问题时，我感到惊讶。 $add = '- 30 days'; echo date('Y-m-01', strtotime($add)); // result is 2017-
我正在使用 javascript 进行练习，我编写了这个脚本来计算 2 个变量的总和，然后在第三个方程中使用这个总和!关于如何完成这项工作的任何想法都将非常有用! First Number:
我有一个来自EAC的提示单和一个包含完整专辑的FLAC文件。 我正在尝试制作一些python脚本来播放文件，因为我需要能够设置在flac文件中开始的位置。 如何从CueSheet格式MM:SS:FF转
这个问题已经有答案了: Adding two numbers concatenates them instead of calculating the sum (24 个回答) 已关闭去年。 我有一个
4000 我需要上面字段 name="quantity" 和 id="price" 中的值，并使用 javascript 函数进行计算，并将其显示在字段 id= 中仅当我单击计算按钮时才显示“总
我是一名优秀的程序员,十分优秀！