src/gmxlib/sfactor.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <ctype.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "string2.h"
  43 #include "futil.h"
  44 #include "maths.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "macros.h"
  48 #include "index.h"
  49 #include "strdb.h"
  50 #include "copyrite.h"
  51 #include "tpxio.h"
  52 #include "typedefs.h"
  53 #include "statutil.h"
  54 #include "oenv.h"
  55 #include "gmxfio.h"
  56 #include "xvgr.h"
  57 #include "matio.h"
  58 #include "gmx_ana.h"
  59 #include "names.h"
  60 #include "sfactor.h"
  61
  62
  63 typedef struct gmx_structurefactors {
  64     int nratoms;
  65     int *p;        /* proton number */
  66     int *n;        /* neutron number */
  67     /* Parameters for the Cromer Mann fit */
  68     real **a;    /* parameter a */
  69     real **b;    /* parameter b */
  70     real *c;       /* parameter c */
  71     char **atomnm;  /* atomname */
  72
  73 } gmx_structurefactors;
  74
  75 typedef struct reduced_atom{
  76     rvec x;
  77     int  t;
  78 } reduced_atom;
  79
  80
  81 typedef struct structure_factor
  82 {
  83   int     n_angles;
  84   int     n_groups;
  85   double  lambda;
  86   double  energy;
  87   double  momentum;
  88   double  ref_k;
  89   double  **F;
  90   int     nSteps;
  91   int     total_n_atoms;
  92 } structure_factor;
  93
  94
  95 extern int * create_indexed_atom_type (reduced_atom_t * atm, int size)
  96 {
  97 /*
  98  * create an index of the atom types found in a  group
  99  * i.e.: for water index_atp[0]=type_number_of_O and
 100  *                 index_atp[1]=type_number_of_H
 101  *
 102  * the last element is set to 0
 103  */
 104     int *index_atp, i, i_tmp, j;
 105
 106     reduced_atom *att=(reduced_atom *)atm;
 107
 108     snew (index_atp, 1);
 109     i_tmp = 1;
 110     index_atp[0] = att[0].t;
 111     for (i = 1; i < size; i++) {
 112         for (j = 0; j < i_tmp; j++)
 113             if (att[i].t == index_atp[j])
 114                 break;
 115         if (j == i_tmp) {       /* i.e. no indexed atom type is  == to atm[i].t */
 116             i_tmp++;
 117             srenew (index_atp, i_tmp * sizeof (int));
 118             index_atp[i_tmp - 1] = att[i].t;
 119         }
 120     }
 121     i_tmp++;
 122     srenew (index_atp, i_tmp * sizeof (int));
 123     index_atp[i_tmp - 1] = 0;
 124     return index_atp;
 125 }
 126
 127
 128
 129 extern t_complex *** rc_tensor_allocation(int x, int y, int z)
 130 {
 131   t_complex ***t;
 132   int i,j;
 133
 134   snew(t,x);
 135   t = (t_complex ***)calloc(x,sizeof(t_complex**));
 136   if(!t) exit(fprintf(stderr,"\nallocation error"));
 137   t[0] = (t_complex **)calloc(x*y,sizeof(t_complex*));
 138   if(!t[0]) exit(fprintf(stderr,"\nallocation error"));
 139   t[0][0] = (t_complex *)calloc(x*y*z,sizeof(t_complex));
 140   if(!t[0][0]) exit(fprintf(stderr,"\nallocation error"));
 141
 142   for( j = 1 ; j < y ; j++)
 143     t[0][j] = t[0][j-1] + z;
 144   for( i = 1 ; i < x ; i++) {
 145     t[i] = t[i-1] + y;
 146     t[i][0] = t[i-1][0] + y*z;
 147     for( j = 1 ; j < y ; j++)
 148       t[i][j] = t[i][j-1] + z;
 149   }
 150   return t;
 151 }
 152
 153
 154 extern void compute_structure_factor (structure_factor_t * sft, matrix box,
 155                                reduced_atom_t * red, int isize, real start_q,
 156                                real end_q, int group,real **sf_table)
 157 {
 158     structure_factor *sf=(structure_factor *)sft;
 159     reduced_atom *redt=(reduced_atom *)red;
 160
 161     t_complex ***tmpSF;
 162     rvec k_factor;
 163     real kdotx, asf, kx, ky, kz, krr;
 164     int kr, maxkx, maxky, maxkz, i, j, k, p, *counter;
 165
 166
 167     k_factor[XX] = 2 * M_PI / box[XX][XX];
 168     k_factor[YY] = 2 * M_PI / box[YY][YY];
 169     k_factor[ZZ] = 2 * M_PI / box[ZZ][ZZ];
 170
 171     maxkx = (int) (end_q / k_factor[XX] + 0.5);
 172     maxky = (int) (end_q / k_factor[YY] + 0.5);
 173     maxkz = (int) (end_q / k_factor[ZZ] + 0.5);
 174
 175     snew (counter, sf->n_angles);
 176
 177     tmpSF = rc_tensor_allocation(maxkx,maxky,maxkz);
 178 /*
 179  * The big loop...
 180  * compute real and imaginary part of the structure factor for every
 181  * (kx,ky,kz))
 182  */
 183     fprintf(stderr,"\n");
 184     for (i = 0; i < maxkx; i++) {
 185         fprintf (stderr,"\rdone %3.1f%%     ", (double)(100.0*(i+1))/maxkx);
 186         kx = i * k_factor[XX];
 187         for (j = 0; j < maxky; j++) {
 188             ky = j * k_factor[YY];
 189             for (k = 0; k < maxkz; k++)
 190                 if (i != 0 || j != 0 || k != 0) {
 191                     kz = k * k_factor[ZZ];
 192                     krr = sqrt (sqr (kx) + sqr (ky) + sqr (kz));
 193                     if (krr >= start_q && krr <= end_q) {
 194                         kr = (int) (krr/sf->ref_k + 0.5);
 195                         if (kr < sf->n_angles) {
 196                             counter[kr]++;  /* will be used for the copmutation
 197                                                of the average*/
 198                             for (p = 0; p < isize; p++) {
 199                                 asf = sf_table[redt[p].t][kr];
 200
 201                                 kdotx = kx * redt[p].x[XX] +
 202                                     ky * redt[p].x[YY] + kz * redt[p].x[ZZ];
 203
 204                                 tmpSF[i][j][k].re += cos (kdotx) * asf;
 205                                 tmpSF[i][j][k].im += sin (kdotx) * asf;
 206                             }
 207                         }
 208                     }
 209                 }
 210         }
 211     }                           /* end loop on i */
 212 /*
 213  *  compute the square modulus of the structure factor, averaging on the surface
 214  *  kx*kx + ky*ky + kz*kz = krr*krr
 215  *  note that this is correct only for a (on the macroscopic scale)
 216  *  isotropic system.
 217  */
 218     for (i = 0; i < maxkx; i++) {
 219         kx = i * k_factor[XX]; for (j = 0; j < maxky; j++) {
 220             ky = j * k_factor[YY]; for (k = 0; k < maxkz; k++) {
 221                 kz = k * k_factor[ZZ]; krr = sqrt (sqr (kx) + sqr (ky)
 222                 + sqr (kz)); if (krr >= start_q && krr <= end_q) {
 223                     kr = (int) (krr / sf->ref_k + 0.5);
 224                         if (kr < sf->n_angles && counter[kr] != 0)
 225                                 sf->F[group][kr] +=
 226                             (sqr (tmpSF[i][j][k].re) +
 227                              sqr (tmpSF[i][j][k].im))/ counter[kr];
 228                 }
 229             }
 230         }
 231     } sfree (counter); free(tmpSF[0][0]); free(tmpSF[0]); free(tmpSF);
 232 }
 233
 234
 235 extern gmx_structurefactors_t *gmx_structurefactors_init(const char *datfn) {
 236
 237         /* Read the database for the structure factor of the different atoms */
 238
 239     FILE *fp;
 240     char line[STRLEN];
 241     gmx_structurefactors *gsf;
 242     double a1,a2,a3,a4,b1,b2,b3,b4,c;
 243     int n,p;
 244     int i;
 245     int nralloc=10;
 246     int line_no;
 247     char atomn[32];
 248     fp=libopen(datfn);
 249     line_no = 0;
 250
 251     snew(gsf,1);
 252
 253     snew(gsf->atomnm,nralloc);
 254     snew(gsf->a,nralloc);
 255     snew(gsf->b,nralloc);
 256     snew(gsf->c,nralloc);
 257     snew(gsf->n,nralloc);
 258     snew(gsf->p,nralloc);
 259     gsf->nratoms=line_no;
 260     while(get_a_line(fp,line,STRLEN)) {
 261         i=line_no;
 262         if (sscanf(line,"%s %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
 263                    atomn,&p,&n,&a1,&a2,&a3,&a4,&b1,&b2,&b3,&b4,&c) == 12) {
 264             gsf->atomnm[i]=strdup(atomn);
 265             gsf->n[i]=n;
 266             gsf->p[i]=p;
 267             snew(gsf->a[i],4);
 268             snew(gsf->b[i],4);
 269             gsf->a[i][0]=a1;
 270             gsf->a[i][1]=a2;
 271             gsf->a[i][2]=a3;
 272             gsf->a[i][3]=a4;
 273             gsf->b[i][0]=b1;
 274             gsf->b[i][1]=b2;
 275             gsf->b[i][2]=b3;
 276             gsf->b[i][3]=b4;
 277             gsf->c[i]=c;
 278             line_no++;
 279             gsf->nratoms=line_no;
 280             if (line_no==nralloc){
 281                 nralloc+=10;
 282                 srenew(gsf->atomnm,nralloc);
 283                 srenew(gsf->a,nralloc);
 284                 srenew(gsf->b,nralloc);
 285                 srenew(gsf->c,nralloc);
 286                 srenew(gsf->n,nralloc);
 287                 srenew(gsf->p,nralloc);
 288             }
 289         }
 290         else
 291             fprintf(stderr,"WARNING: Error in file %s at line %d ignored\n",
 292                     datfn,line_no);
 293     }
 294
 295     srenew(gsf->atomnm,gsf->nratoms);
 296     srenew(gsf->a,gsf->nratoms);
 297     srenew(gsf->b,gsf->nratoms);
 298     srenew(gsf->c,gsf->nratoms);
 299     srenew(gsf->n,gsf->nratoms);
 300     srenew(gsf->p,gsf->nratoms);
 301
 302     fclose(fp);
 303
 304     return (gmx_structurefactors_t *) gsf;
 305
 306 }
 307
 308
 309 extern void rearrange_atoms (reduced_atom_t * positions, t_trxframe *fr, atom_id * index,
 310                       int isize, t_topology * top, bool flag,gmx_structurefactors_t *gsf)
 311 /* given the group's index, return the (continuous) array of atoms */
 312 {
 313   int i;
 314
 315   reduced_atom *pos=(reduced_atom *)positions;
 316
 317   if (flag)
 318     for (i = 0; i < isize; i++)
 319       pos[i].t =
 320         return_atom_type (*(top->atoms.atomname[index[i]]),gsf);
 321   for (i = 0; i < isize; i++)
 322     copy_rvec (fr->x[index[i]], pos[i].x);
 323
 324    positions=(reduced_atom_t *)pos;
 325 }
 326
 327
 328 extern int return_atom_type (const char *name,gmx_structurefactors_t *gsf)
 329 {
 330   typedef struct {
 331     const char *name;
 332     int  nh;
 333   } t_united_h;
 334   t_united_h uh[] = {
 335     { "CH1", 1 }, { "CH2", 2 }, { "CH3", 3 },
 336     { "CS1", 1 }, { "CS2", 2 }, { "CS3", 3 },
 337     { "CP1", 1 }, { "CP2", 2 }, { "CP3", 3 }
 338   };
 339   int i,cnt=0;
 340   int *tndx;
 341   int nrc;
 342   int fndx=0;
 343   int NCMT;
 344
 345   gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
 346
 347   NCMT=gsft->nratoms;
 348
 349   snew(tndx,NCMT);
 350
 351   for(i=0; (i<asize(uh)); i++)
 352     if (strcmp(name,uh[i].name) == 0)
 353       return NCMT-1+uh[i].nh;
 354
 355   for(i=0; (i<NCMT); i++){
 356     if (strncmp (name, gsft->atomnm[i],strlen(gsft->atomnm[i])) == 0){
 357       tndx[cnt]=i;
 358       cnt++;
 359     }
 360   }
 361
 362   if (cnt==0)
 363   gmx_fatal(FARGS,"\nError: atom (%s) not in list (%d types checked)!\n",
 364             name,i);
 365   else{
 366           nrc=0;
 367           for(i=0;i<cnt;i++){
 368                   if(strlen(gsft->atomnm[tndx[i]])>(size_t)nrc){
 369                           nrc=strlen(gsft->atomnm[tndx[i]]);
 370                       fndx=tndx[i];
 371                   }
 372           }
 373
 374           return fndx;
 375   }
 376
 377   return 0;
 378 }
 379
 380 extern int gmx_structurefactors_get_sf(gmx_structurefactors_t *gsf, int elem, real a[4], real b[4], real *c){
 381
 382         int success;
 383         int i;
 384     gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
 385         success=0;
 386
 387         for(i=0;i<4;i++){
 388             a[i]=gsft->a[elem][i];
 389             b[i]=gsft->b[elem][i];
 390             *c=gsft->c[elem];
 391         }
 392
 393         success+=1;
 394         return success;
 395 }
 396
 397 int atp_size (int *index_atp)
 398 {
 399     int i = 0;
 400
 401     while (index_atp[i])
 402         i++;
 403     return i;
 404 }
 405
 406
 407 extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
 408                              const char* fnXVG, const char *fnTRX,
 409                              const char* fnDAT,
 410                              real start_q,real end_q,
 411                              real energy,int ng,const output_env_t oenv)
 412 {
 413     int i,*isize,flags = TRX_READ_X,**index_atp;
 414     t_trxstatus *status;
 415     char **grpname,title[STRLEN];
 416     atom_id **index;
 417     t_topology top;
 418     int ePBC;
 419     t_trxframe fr;
 420     reduced_atom_t **red;
 421     structure_factor *sf;
 422     rvec *xtop;
 423     real **sf_table;
 424     int nsftable;
 425     matrix box;
 426     double r_tmp;
 427
 428     gmx_structurefactors_t *gmx_sf;
 429     real *a,*b,c;
 430     int success;
 431
 432     snew(a,4);
 433     snew(b,4);
 434
 435
 436     gmx_sf=gmx_structurefactors_init(fnDAT);
 437
 438     success=gmx_structurefactors_get_sf(gmx_sf,0, a, b, &c);
 439
 440     snew (sf, 1);
 441     sf->energy = energy;
 442
 443     /* Read the topology informations */
 444     read_tps_conf (fnTPS, title, &top, &ePBC, &xtop, NULL, box, TRUE);
 445     sfree (xtop);
 446
 447     /* groups stuff... */
 448     snew (isize, ng);
 449     snew (index, ng);
 450     snew (grpname, ng);
 451
 452     fprintf (stderr, "\nSelect %d group%s\n", ng,
 453              ng == 1 ? "" : "s");
 454     if (fnTPS)
 455         get_index (&top.atoms, fnNDX, ng, isize, index, grpname);
 456     else
 457         rd_index (fnNDX, ng, isize, index, grpname);
 458
 459     /* The first time we read data is a little special */
 460     read_first_frame (oenv,&status, fnTRX, &fr, flags);
 461
 462     sf->total_n_atoms = fr.natoms;
 463
 464     snew (red, ng);
 465     snew (index_atp, ng);
 466
 467     r_tmp = max (box[XX][XX], box[YY][YY]);
 468     r_tmp = (double) max (box[ZZ][ZZ], r_tmp);
 469
 470     sf->ref_k = (2.0 * M_PI) / (r_tmp);
 471     /* ref_k will be the reference momentum unit */
 472     sf->n_angles = (int) (end_q / sf->ref_k + 0.5);
 473
 474     snew (sf->F, ng);
 475     for (i = 0; i < ng; i++)
 476         snew (sf->F[i], sf->n_angles);
 477     for (i = 0; i < ng; i++) {
 478         snew (red[i], isize[i]);
 479         rearrange_atoms (red[i], &fr, index[i], isize[i], &top, TRUE,gmx_sf);
 480         index_atp[i] = create_indexed_atom_type (red[i], isize[i]);
 481     }
 482
 483     sf_table = compute_scattering_factor_table (gmx_sf,(structure_factor_t *)sf,&nsftable);
 484
 485
 486     /* This is the main loop over frames */
 487
 488     do {
 489         sf->nSteps++;
 490         for (i = 0; i < ng; i++) {
 491             rearrange_atoms (red[i], &fr, index[i], isize[i], &top,FALSE,gmx_sf);
 492
 493             compute_structure_factor ((structure_factor_t *)sf, box, red[i], isize[i],
 494                                       start_q, end_q, i, sf_table);
 495         }
 496     }
 497
 498     while (read_next_frame (oenv,status, &fr));
 499
 500     save_data ((structure_factor_t *)sf, fnXVG, ng, start_q, end_q,oenv);
 501
 502
 503     sfree(a);
 504     sfree(b);
 505
 506     gmx_structurefactors_done(gmx_sf);
 507
 508     return 0;
 509 }
 510
 511
 512 extern void save_data (structure_factor_t *sft, const char *file, int ngrps,
 513                 real start_q, real end_q, const output_env_t oenv)
 514 {
 515
 516     FILE *fp;
 517     int i, g = 0;
 518     double *tmp, polarization_factor, A;
 519
 520     structure_factor *sf=(structure_factor *)sft;
 521
 522     fp = xvgropen (file, "Scattering Intensity", "q (1/nm)",
 523                    "Intensity (a.u.)",oenv);
 524
 525     snew (tmp, ngrps);
 526
 527     for (g = 0; g < ngrps; g++)
 528         for (i = 0; i < sf->n_angles; i++) {
 529
 530 /*
 531  *          theta is half the angle between incoming and scattered vectors.
 532  *
 533  *          polar. fact. = 0.5*(1+cos^2(2*theta)) = 1 - 0.5 * sin^2(2*theta)
 534  *
 535  *          sin(theta) = q/(2k) := A  ->  sin^2(theta) = 4*A^2 (1-A^2) ->
 536  *          -> 0.5*(1+cos^2(2*theta)) = 1 - 2 A^2 (1-A^2)
 537  */
 538             A = (double) (i * sf->ref_k) / (2.0 * sf->momentum);
 539             polarization_factor = 1 - 2.0 * sqr (A) * (1 - sqr (A));
 540             sf->F[g][i] *= polarization_factor;
 541         }
 542     for (i = 0; i < sf->n_angles; i++) {
 543         if (i * sf->ref_k >= start_q && i * sf->ref_k <= end_q) {
 544             fprintf (fp, "%10.5f  ", i * sf->ref_k);
 545             for (g = 0; g < ngrps; g++)
 546                fprintf (fp, "  %10.5f ", (sf->F[g][i]) /( sf->total_n_atoms*
 547                                                           sf->nSteps));
 548             fprintf (fp, "\n");
 549         }
 550     }
 551
 552     ffclose (fp);
 553 }
 554
 555
 556 extern double CMSF (gmx_structurefactors_t *gsf,int type,int nh,double lambda, double sin_theta)
 557 /*
 558  * return Cromer-Mann fit for the atomic scattering factor:
 559  * sin_theta is the sine of half the angle between incoming and scattered
 560  * vectors. See g_sq.h for a short description of CM fit.
 561  */
 562 {
 563   int i,success;
 564   double tmp = 0.0, k2;
 565   real *a,*b;
 566   real c;
 567
 568   snew(a,4);
 569   snew(b,4);
 570
 571  /*
 572   *
 573   * f0[k] = c + [SUM a_i*EXP(-b_i*(k^2)) ]
 574   *             i=1,4
 575   */
 576
 577   /*
 578    *  united atoms case
 579    *  CH2 / CH3 groups
 580    */
 581   if (nh > 0) {
 582     tmp = (CMSF (gsf,return_atom_type ("C",gsf),0,lambda, sin_theta) +
 583            nh*CMSF (gsf,return_atom_type ("H",gsf),0,lambda, sin_theta));
 584   }
 585   /* all atom case */
 586   else {
 587     k2 = (sqr (sin_theta) / sqr (10.0 * lambda));
 588     success=gmx_structurefactors_get_sf(gsf,type,a,b,&c);
 589     tmp = c;
 590     for (i = 0; (i < 4); i++)
 591       tmp += a[i] * exp (-b[i] * k2);
 592   }
 593   return tmp;
 594 }
 595
 596
 597
 598 extern real **gmx_structurefactors_table(gmx_structurefactors_t *gsf,real momentum, real ref_k, real lambda, int n_angles){
 599
 600         int NCMT;
 601         int nsftable;
 602         int i,j;
 603         double q,sin_theta;
 604     real **sf_table;
 605     gmx_structurefactors *gsft=(gmx_structurefactors *)gsf;
 606
 607     NCMT=gsft->nratoms;
 608         nsftable = NCMT+3;
 609
 610     snew (sf_table,nsftable);
 611     for (i = 0; (i < nsftable); i++) {
 612         snew (sf_table[i], n_angles);
 613         for (j = 0; j < n_angles; j++) {
 614                 q = ((double) j * ref_k);
 615                 /* theta is half the angle between incoming
 616                            and scattered wavevectors. */
 617                 sin_theta = q / (2.0 * momentum);
 618                 if (i < NCMT){
 619                     sf_table[i][j] = CMSF (gsf,i,0,lambda, sin_theta);
 620                 }
 621                 else
 622                         sf_table[i][j] = CMSF (gsf,i,i-NCMT+1,lambda, sin_theta);
 623         }
 624     }
 625     return sf_table;
 626 }
 627
 628 extern void gmx_structurefactors_done(gmx_structurefactors_t *gsf){
 629
 630         int i;
 631         gmx_structurefactors *sf;
 632         sf=(gmx_structurefactors *) gsf;
 633
 634         for(i=0;i<sf->nratoms;i++){
 635         sfree(sf->a[i]);
 636                 sfree(sf->b[i]);
 637                 sfree(sf->atomnm[i]);
 638         }
 639
 640         sfree(sf->a);
 641         sfree(sf->b);
 642         sfree(sf->atomnm);
 643         sfree(sf->p);
 644         sfree(sf->n);
 645         sfree(sf->c);
 646
 647         sfree(sf);
 648
 649 }
 650
 651 extern real **compute_scattering_factor_table (gmx_structurefactors_t *gsf,structure_factor_t *sft,int *nsftable)
 652 {
 653 /*
 654  *  this function build up a table of scattering factors for every atom
 655  *  type and for every scattering angle.
 656  */
 657
 658     double hc=1239.842;
 659     real ** sf_table;
 660
 661     structure_factor *sf=(structure_factor *)sft;
 662
 663
 664     /* \hbar \omega \lambda = hc = 1239.842 eV * nm */
 665     sf->momentum = ((double) (2. * 1000.0 * M_PI * sf->energy) / hc);
 666     sf->lambda = hc / (1000.0 * sf->energy);
 667     fprintf (stderr, "\nwavelenght = %f nm\n", sf->lambda);
 668
 669     sf_table=gmx_structurefactors_table(gsf,sf->momentum,sf->ref_k,sf->lambda,sf->n_angles);
 670
 671     return sf_table;
 672 }
 673
 674