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