1/*******************************************************************************
2* McStas, neutron ray-tracing package
3*         Copyright 1997-2003, All rights reserved
4*         Risoe National Laboratory, Roskilde, Denmark
5*         Institut Laue Langevin, Grenoble, France
6*
7* Component: SANSPDBFast
8*
9* %I
10* Written by: Martin Cramer Pedersen (mcpe@nbi.dk) and Søren Kynde (kynde@nbi.dk)
11* Date: October 17, 2012
12* Origin: KU-Science
13*
14* A sample describing a thin solution of proteins using linear interpolation
15* to increase computational speed. This components must be compiled with the
16* -lgsl and -lgslcblas flags (and possibly linked to the appropriate libraries).
17*
18* %D
19* This components expands the formfactor amplitude of the protein on spherical
20* harmonics and computes the scattering profile using these. The expansion is
21* done on amino-acid level and does not take hydration layer into account.
22* The component must have a valid .pdb-file as an argument.
23*
24* %P
25* RhoSolvent: [AA]               Scattering length density of the buffer - default is 100% D2O.
26* Concentration: [mM]            Concentration of sample.
27* AbsorptionCrosssection: [1/m]  Absorption cross section of the sample.
28* xwidth: [m]                    Dimension of component in the x-direction.
29* yheight: [m]                   Dimension of component in the y-direction.
30* zdepth: [m]                    Dimension of component in the z-direction.
31* SampleToDetectorDistance: [m]  Distance from sample to detector (for focusing the scattered neutrons).
32* DetectorRadius: [m]            Radius of the detector (for focusing the scattered neutrons).
33* qMin: [AA^-1]                  Lowest q-value, for which a point is generated in the scattering profile
34* qMax: [AA^-1]                  Highest q-value, for which a point is generated in the scattering profile
35* NumberOfQBins: []		Number of points generated in initalscattering profile.
36* PDBFilepath: []		Path to the file describing the high resolution structure of the protein.
37*
38* %E
39*******************************************************************************/
40
41DEFINE COMPONENT SANSPDBFast
42
43DEFINITION PARAMETERS (NumberOfQBins = 200)
44
45SETTING PARAMETERS (RhoSolvent = 6.4e-14, Concentration = 0.01, AbsorptionCrosssection = 0.0,
46xwidth, yheight, zdepth,
47SampleToDetectorDistance, DetectorRadius,
48qMin = 0.001, qMax = 0.5,
49string PDBFilepath = "PDBfile.pdb")
50
51OUTPUT PARAMETERS ()
52
53DECLARE
54%{
55// Declarations
56double Absorption;
57double q;
58double NumberDensity;
59const int OrderOfHarmonics;
60
61// Arrays for storing q and I(q)
62	double * qArray;
63	double * IArray;
64
65	// GSL libraries
66	#include <gsl/gsl_sf_legendre.h>
67	#include <gsl/gsl_sf_bessel.h>
68	#include <complex.h>
69%}
70
71INITIALIZE
72%{
73	OrderOfHarmonics = 21;
74
75	// Rescale concentration into number of aggregates per m^3 times 10^-4
76	NumberDensity = Concentration * 6.02214129e19;
77
78	// Structs describing the protein
79	struct Bead
80	{
81		double x;
82		double y;
83		double z;
84
85		double xv;
86		double yv;
87		double zv;
88
89		double xa;
90		double ya;
91		double za;
92
93		char *NameOfResidue;
94
95		double Volume;
96		double ScatteringLength;
97	};
98	typedef struct Bead BeadStruct;
99
100	struct Protein
101	{
102		BeadStruct *Beads;
103		int NumberOfResidues;
104	};
105	typedef struct Protein ProteinStruct;
106
107	// Simple mathematical functions
108	int Sign(double x) {
109		int Sign;
110
111		if (x > 0) {
112			Sign = 1;
113		} else if (x < 0) {
114			Sign = -1;
115		} else {
116			Sign = 0;
117		}
118
119		return Sign;
120	}
121
122	double complex Polar(double R, double Concentration) {
123		double complex Polar;
124
125		Polar = R * (cos(Concentration) + _Complex_I * sin(Concentration));
126
127		return Polar;
128	}
129
130	// Function used to determine the number of residues in the .pdb-file
131	int CountResidues(char PDBFilepath[256])
132	{
133		// Declarations
134		double Dummy1;
135		double Dummy2;
136		double Dummy3;
137		char Line[256];
138		char DummyChar;
139		char Atom;
140		int NumberOfResidues = 0;
141		int ResidueID;
142		int PreviousResidueID = 0;
143		FILE *PDBFile;
144
145		// I/O
146		PDBFile = fopen(PDBFilepath, "r");
147
148		if (PDBFile == NULL) {
149			printf("Cannot open %s... \n", PDBFilepath);
150		} else {
151			printf(".pdb-file located... \n");
152		}
153
154		while (fgets(Line, sizeof(Line), PDBFile) != NULL) {
155		    ResidueID = 0;
156
157		    if (sscanf(Line, "ATOM%*18c%d%*4c%lf%lf%lf", &ResidueID, &Dummy1, &Dummy2, &Dummy3) == 4) {
158
159		        if (ResidueID != PreviousResidueID && ResidueID != 0) {
160		            ++NumberOfResidues;
161				}
162
163		        PreviousResidueID = ResidueID;
164		    }
165		}
166
167		fclose(PDBFile);
168		printf("Calculating scattering from %d residues...\n", NumberOfResidues);
169
170        return NumberOfResidues;
171	}
172
173	// Initialization of the arrays
174	double complex ** ComplexMatrix(int A, int B)
175	{
176		int i;
177		int j;
178		double complex ** Matrix;
179
180		Matrix = (double complex **) malloc(A * sizeof(double complex *));
181
182		for (i = 0; i < A; ++i) {
183		    Matrix[i] = (double complex *) malloc(B * sizeof(double complex));
184
185		    for (j = 0; j < B; ++j) {
186		        Matrix[i][j] = 0.0;
187		    }
188		}
189
190		return Matrix;
191	}
192
193	void InitializeProtein(ProteinStruct *Protein, int NumberOfResiduesInFile)
194	{
195		Protein->NumberOfResidues = NumberOfResiduesInFile;
196		Protein->Beads = calloc((size_t) NumberOfResiduesInFile, sizeof(BeadStruct));
197	}
198
199	void InitializeArray(int Size, double** Array)
200	{
201		// Declarations
202		int i;
203		double * DummyArray;
204
205		// Initialization
206		DummyArray = (double *) calloc(Size, sizeof(double));
207
208		for (i = 0; i < Size; ++i) {
209			DummyArray[i] = 0.0;
210		}
211
212        *Array = DummyArray;
213	}
214
215	// Function used to read .pdb-file
216	void ReadAminoPDB(char PDBFilename[256], ProteinStruct *Protein)
217	{
218		// Declarations and input
219		int NumberOfResidues = Protein->NumberOfResidues;
220		BeadStruct *Residue = Protein->Beads;
221		FILE *PDBFile;
222
223		int i = 0;
224		int PreviousResidueID = 0;
225		int ResidueID = 0;
226
227		double Weight = 0.0;
228		double W = 0.0;
229
230		double Aweight = 0.0;
231		double A = 0.0;
232
233		double x;
234		double y;
235		double z;
236
237		double X = 0.0;
238		double Y = 0.0;
239		double Z = 0.0;
240
241		double XA = 0.0;
242		double YA = 0.0;
243		double ZA = 0.0;
244
245		char Atom;
246
247		char Line[256];
248		char Buffer[256];
249		char DummyChar;
250
251		// Atomic weighing factors
252		const double WH = 5.15;
253		const double WC = 16.44;
254		const double WN = 2.49;
255		const double WO = 9.13;
256		const double WS = 19.86;
257		const double WP = 5.73;
258
259		// Scattering lengths
260		const double AH =  - 3.741e-15;
261		const double AD =    6.674e-15;
262		const double AC =    6.648e-15;
263		const double AN =    9.360e-15;
264		const double AO =    5.805e-15;
265		const double AP =    5.130e-15;
266		const double AS =    2.847e-15;
267
268		// Program
269		if ((PDBFile = fopen(PDBFilename, "r")) == 0) {
270		    printf("Cannot open file: %s. \n", PDBFilename);
271		}
272
273		Residue[i].NameOfResidue = (char *) calloc(3, sizeof(char));
274
275		while (fgets(Buffer, sizeof(Buffer), PDBFile) != NULL) {
276		    Atom = 0;
277			ResidueID = 0;
278
279			sscanf(Buffer,"ATOM%*9c%c%*8c%d%*4c%lf%lf%lf%*23c%c", &DummyChar, &ResidueID, &x, &y, &z, &Atom);
280
281		    if (ResidueID != PreviousResidueID && ResidueID != 0) {
282
283		        if (PreviousResidueID != 0) {
284
285					// Assign center of scattering
286				    Residue[i].xv = X / Weight;
287					Residue[i].yv = Y / Weight;
288					Residue[i].zv = Z / Weight;
289
290					// Assign center of mass
291				    Residue[i].x = XA / Aweight;
292					Residue[i].y = YA / Aweight;
293					Residue[i].z = ZA / Aweight;
294
295					// Other residue attributes
296				    Residue[i].Volume = Weight;
297				    Residue[i].ScatteringLength = Aweight;
298
299				    X = 0.0;
300					Y = 0.0;
301					Z = 0.0;
302					Weight = 0.0;
303
304					XA = 0.0;
305					YA = 0.0;
306					ZA = 0.0;
307					Aweight = 0.0;
308
309				    ++i;
310
311				    if (i < NumberOfResidues) {
312				        Residue[i].NameOfResidue = (char *) calloc(3, sizeof(char));
313				    }
314
315				}
316
317		        PreviousResidueID = ResidueID;
318		    }
319
320			// Finish the final amino acid
321			if (i == NumberOfResidues - 1) {
322				Residue[i].xv = X / Weight;
323				Residue[i].yv = Y / Weight;
324				Residue[i].zv = Z / Weight;
325
326				// Assign center of mass
327				Residue[i].x = XA / Aweight;
328				Residue[i].y = YA / Aweight;
329				Residue[i].z = ZA / Aweight;
330
331				// Other residue attributes
332				Residue[i].Volume = Weight;
333				Residue[i].ScatteringLength = Aweight;
334			}
335
336		    sscanf(Buffer, "ATOM%*9cCA%*2c%*s%*10c%lf%lf%lf", &Residue[i].xa, &Residue[i].ya, &Residue[i].za);
337		    sscanf(Buffer, "ATOM%*13c%s", Residue[i].NameOfResidue);
338
339		    switch(Atom) {
340		        case 'C':
341		            A = AC;
342		            W = WC;
343		        break;
344
345		        case 'N':
346		            A = AN;
347		            W = WN;
348		        break;
349
350		        case 'O':
351		            A = AO;
352		            W = WO;
353		        break;
354
355		        case 'S':
356		            A = AS;
357		            W = WS;
358		        break;
359
360		        case 'H':
361		            A = AH;
362		            W = WH;
363		        break;
364
365		        case 'P':
366		            A = AP;
367		            W = WP;
368		        break;
369
370		        default:
371		            A = 0.0;
372		            W = 0.0;
373		    }
374
375		    Weight += W;
376		    Aweight += A;
377
378		    X += W * x;
379		    Y += W * y;
380		    Z += W * z;
381
382		    XA += A * x;
383		    YA += A * y;
384		    ZA += A * z;
385		}
386
387		fclose(PDBFile);
388	}
389
390	// Function used to expand protein structure on spherical harmonics
391	void ExpandStructure(double complex ** Matrix, ProteinStruct *Protein, int ResidueID, double q, double RhoSolvent)
392	{
393		// Dummy integers
394		int i;
395		int j;
396
397		// Arrays used for storing Legendre and Bessel function values
398		double Legendre[OrderOfHarmonics + 1];
399		double Bessel[OrderOfHarmonics + 1];
400
401		// Residue information
402		const double Volume = Protein->Beads[ResidueID].Volume;
403		const double DeltaRhoProtein = Protein->Beads[ResidueID].ScatteringLength - Volume * RhoSolvent;
404
405		const double x = (Protein->Beads[ResidueID].x * Protein->Beads[ResidueID].ScatteringLength -
406						  RhoSolvent * Volume * Protein->Beads[ResidueID].xv) / DeltaRhoProtein;
407
408		const double y = (Protein->Beads[ResidueID].y * Protein->Beads[ResidueID].ScatteringLength -
409						  RhoSolvent * Volume * Protein->Beads[ResidueID].yv) / DeltaRhoProtein;
410
411		const double z = (Protein->Beads[ResidueID].z * Protein->Beads[ResidueID].ScatteringLength -
412						  RhoSolvent * Volume * Protein->Beads[ResidueID].zv) / DeltaRhoProtein;
413
414		// Convert bead position to spherical coordinates
415		const double Radius = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
416		const double Theta  = acos(z / Radius);
417		const double Concentration    = acos(x / (Radius * sin(Theta))) * Sign(y);
418
419		// Expand protein structure on harmonics
420		gsl_sf_bessel_jl_array(OrderOfHarmonics, q * Radius, Bessel);
421
422		for (i = 0; i <= OrderOfHarmonics; ++i) {
423		    gsl_sf_legendre_sphPlm_array(OrderOfHarmonics, i, cos(Theta), &Legendre[i]);
424
425		    for(j = i; j <= OrderOfHarmonics; ++j) {
426		        Matrix[j][i] += sqrt(4.0 * PI) * cpow(_Complex_I, j) * DeltaRhoProtein * Bessel[j] * Legendre[j] * Polar(1.0, -i * Concentration);
427		    }
428		}
429	}
430
431	// Function used to generate intensity for a given q-value
432	double ComputeIntensity(double complex ** Matrix, int Size)
433	{
434		// Declarations
435		int i;
436		int j;
437		double Intensity = 0.0;
438
439		// Computation
440		for (i = 0; i <= Size; ++i) {
441
442		    for (j = 0; j <= i; ++j) {
443		        Intensity += ((j > 0) + 1.0) * pow(cabs(Matrix[i][j]), 2);
444		    }
445		}
446
447		return Intensity;
448	}
449
450	// Function used to reinitialize a matrix
451	void ResetMatrix(double complex ** Matrix, int Size)
452	{
453		int i;
454		int j;
455
456		for (i = 0; i <= Size; ++i) {
457
458		    for(j = 0; j <= Size; ++j) {
459		        Matrix[i][j] = 0.0;
460		    }
461		}
462	}
463
464	// Standard sample handling
465	if (!xwidth || !yheight || !zdepth) {
466		printf("%s: Sample has no volume - check parameters...! \n", NAME_CURRENT_COMP);
467	}
468
469	Absorption = AbsorptionCrosssection;
470
471	// Declarations
472	int i;
473	int j;
474	double qDummy;
475	const double qStep = (qMax - qMin) / (1.0 * NumberOfQBins);
476
477	// Initialize matrix
478	double complex ** Matrix = ComplexMatrix(OrderOfHarmonics + 1, OrderOfHarmonics + 1);
479
480	// Initialize protein
481	int NumberOfResidues;
482	ProteinStruct Protein;
483
484	printf("Initializing protein structure...\n");
485    NumberOfResidues = CountResidues(PDBFilepath);
486	InitializeProtein(&Protein, NumberOfResidues);
487
488	printf("Initializing arrays...\n");
489	InitializeArray(NumberOfQBins, &qArray);
490	InitializeArray(NumberOfQBins, &IArray);
491
492	printf("Creating protein structure...\n");
493	ReadAminoPDB(PDBFilepath, &Protein);
494
495	// Computing scattering profile
496	printf("Computing scattering from protein...\n");
497	for (i = 0; i < NumberOfQBins; ++i) {
498		ResetMatrix(Matrix, OrderOfHarmonics);
499
500		qDummy = qMin + (i + 0.5) * qStep;
501
502		for (j = 0; j < NumberOfResidues; ++j) {
503			ExpandStructure(Matrix, &Protein, j, qDummy, RhoSolvent);
504		}
505
506		qArray[i] = qDummy;
507	 	IArray[i] = ComputeIntensity(Matrix, OrderOfHarmonics);
508	}
509
510	printf("Initializations complete...\n");
511%}
512
513TRACE
514%{
515	// Declarations
516	double t0;
517	double t1;
518	double l_full;
519	double l;
520	double l1;
521	double Intensity;
522	double Weight;
523	double IntensityPart;
524	double SolidAngle;
525	double qx;
526	double qy;
527	double qz;
528	double v;
529	double dt;
530	double vx_i;
531	double vy_i;
532	double vz_i;
533	char Intersect = 0;
534	double Slope;
535	double Offset;
536	int i;
537
538	// Computation
539	Intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
540
541	if (Intersect) {
542
543		if (t0 < 0.0) {
544			fprintf(stderr, "Neutron already inside sample %s - absorbing...\n", NAME_CURRENT_COMP);
545			ABSORB;
546    	}
547
548		// Compute properties of neutron
549		v = sqrt(pow(vx, 2) + pow(vy, 2) + pow(vz, 2));
550		l_full = v * (t1 - t0);
551		dt = rand01() * (t1 - t0) + t0;
552		PROP_DT(dt);
553	    l = v * (dt - t0);
554
555		// Store properties of incoming neutron
556		vx_i = vx;
557		vy_i = vy;
558		vz_i = vz;
559
560		// Generate new direction of neutron
561		randvec_target_circle(&vx, &vy, &vz, &SolidAngle, 0, 0, SampleToDetectorDistance, DetectorRadius);
562
563		NORM(vx, vy, vz);
564
565		vx *= v;
566		vy *= v;
567		vz *= v;
568
569		// Compute q
570		qx = V2K * (vx_i - vx);
571		qy = V2K * (vy_i - vy);
572		qz = V2K * (vz_i - vz);
573
574		q = sqrt(pow(qx, 2) + pow(qy, 2) + pow(qz, 2));
575
576		// Discard neutron, if q is out of range
577		if ((q < qArray[0]) || (q > qArray[NumberOfQBins - 1])) {
578		    ABSORB;
579		}
580
581		// Find the first value of q in the curve larger than that of the neutron
582		i = 1;
583
584		while (q > qArray[i]) {
585			++i;
586		}
587
588		// Do a linear interpolation
589		l1 = v * t1;
590
591		Slope = (IArray[i] - IArray[i - 1]) / (qArray[i] - qArray[i - 1]);
592		Offset = IArray[i] - Slope * qArray[i];
593
594		Intensity = (Slope * q + Offset);
595
596		p *= l_full * SolidAngle / (4.0 * PI) * NumberDensity * Intensity * exp(- Absorption * (l + l1) / v);
597
598		SCATTER;
599	}
600%}
601
602MCDISPLAY
603%{
604
605	box(0, 0, 0, xwidth, yheight, zdepth);
606%}
607
608END
609