src/gromacs/gmxana/sfactor.cpp

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