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