src/gromacs/gmxana/sfactor.c

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