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