-
Notifications
You must be signed in to change notification settings - Fork 0
/
tpfft.c
109 lines (93 loc) · 2.59 KB
/
tpfft.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include <stdlib.h>
#include <string.h>
#include "dcd.h"
#include "psf.h"
#include "betaz.h"
#include "pfft.h"
#include "readArray.h"
void dumpRes(struct result res) {
printf("q,I(q)\n");
for (int i=0; i<res.nqs; i++)
printf("%lf,%lf\n",res.qs[i],res.is[i]);
}
void rewrap(float * xs, int n, double l) {
for (int i=0; i<n; i++) {
while(xs[i] >= l)
xs[i]-=l;
while(xs[i] < 0)
xs[i]+=l;
}
}
int main(int argc, char *argv[]) {
//get scattering lengths
fprintf(stderr,"Reading structure file.\n");
struct psf p = readPSF(argv[1]);
int natom=p.natom;
double bs[natom];
for(int i=0; i<natom; i++)
bs[i]=p.atoms[i].b;
freePSF(p);
//prep coords arrays
fprintf(stderr,"Opening DCD.\n");
struct dcd * d = openDCD(argv[2]);
double uc[6];
float xs[natom];
float ys[natom];
float zs[natom];
struct betaz * bz = NULL;
//build beta(z)
fprintf(stderr,"Generating beta(z)\n");
for (int i=0; i<getNFrames(d); i++) {
getUnitCell(d,uc);
getCoords(d,xs,ys,zs);
bz = runAvgBetaZ(bs,zs,natom, uc[0]*uc[1], uc[2], 0.5, bz);
nextFrame(d);
}
//get solvent NSLD from beta(z)
double bw=0;
if(bz->maxbin>30)
bw = (bz->beta[-bz->maxbin] +
bz->beta[-bz->maxbin + 1] +
bz->beta[-bz->maxbin + 2] +
bz->beta[-bz->maxbin + 3] +
bz->beta[ bz->maxbin - 3] +
bz->beta[ bz->maxbin - 2] +
bz->beta[ bz->maxbin - 1] +
bz->beta[ bz->maxbin])/8.0;
else
bw = (bz->beta[-bz->maxbin] + bz->beta[ bz->maxbin])/2.0;
goToFrame(d,0);
//get qs
fprintf(stderr,"Reading qs\n");
double * qs = malloc(1024*1024*sizeof(double));
int nqs = readDoubles(argv[3],qs);
while(qs[nqs-1]==0)
nqs--;
qs = realloc(qs,nqs*sizeof(double));
double * is = calloc(nqs,sizeof(double));
double * tis = calloc(nqs,sizeof(double));
struct partsys ps = {.natom = natom, .xs = xs, .ys = ys, .zs = zs, .bs = bs, .uc = uc};
struct result res = {.nqs = nqs, .qs = qs, .is = tis};
double ds = 1; //1 angstom bin size
int rci = 15; //30 angstom cutoff
int simax = 199; //200 angstom limit
fprintf(stderr,"Running pfft.\n");
for (int i=0; i<getNFrames(d); i++) {
getUnitCell(d,uc);
getCoords(d,xs,ys,zs);
rewrap(xs,natom,uc[0]);
rewrap(ys,natom,uc[1]);
rewrap(zs,natom,uc[2]);
memset(res.is,0,res.nqs*sizeof(double));
pfft(ps,bz,bw,rci,ds,simax,res);
for(int qi=0; qi<res.nqs; qi++) {
is[qi]+= (res.is[qi]-is[qi])/(i + 1);
}
}
free(tis);
dumpRes(res);
closeDCD(d);
free(qs);
free(is);
freeBetaZ(bz);
}