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