src/gromacs/gmxana/nsfactor.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <string.h>
  40
  41 #include "gromacs/random/random.h"
  42 #include "gromacs/utility/smalloc.h"
  43 #include "sysstuff.h"
  44 #include "vec.h"
  45 #include "nsfactor.h"
  46
  47 #include "gromacs/fileio/futil.h"
  48 #include "gromacs/fileio/strdb.h"
  49 #include "gmx_fatal.h"
  50 #include "gromacs/utility/gmxomp.h"
  51
  52 void check_binwidth(real binwidth)
  53 {
  54     real smallest_bin = 0.1;
  55     if (binwidth < smallest_bin)
  56     {
  57         gmx_fatal(FARGS, "Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
  58     }
  59 }
  60
  61 void check_mcover(real mcover)
  62 {
  63     if (mcover > 1.0)
  64     {
  65         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  66     }
  67     else if ((mcover < 0)&(mcover != -1))
  68     {
  69         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  70     }
  71     else
  72     {
  73         return;
  74     }
  75 }
  76
  77 void normalize_probability(int n, double *a)
  78 {
  79     int    i;
  80     double norm = 0.0;
  81     for (i = 0; i < n; i++)
  82     {
  83         norm += a[i];
  84     }
  85     for (i = 0; i < n; i++)
  86     {
  87         a[i] /= norm;
  88     }
  89 }
  90
  91 gmx_neutron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn)
  92 {
  93     /* read nsfactor.dat */
  94     FILE    *fp;
  95     char     line[STRLEN];
  96     int      nralloc = 10;
  97     int      n, p;
  98     int      i, line_no;
  99     char     atomnm[8];
 100     double   slength;
 101     gmx_neutron_atomic_structurefactors_t   *gnsf;
 102
 103     fp      = libopen(datfn);
 104     line_no = 0;
 105     /* allocate memory for structure */
 106     snew(gnsf, nralloc);
 107     snew(gnsf->atomnm, nralloc);
 108     snew(gnsf->p, nralloc);
 109     snew(gnsf->n, nralloc);
 110     snew(gnsf->slength, nralloc);
 111
 112     gnsf->nratoms = line_no;
 113
 114     while (get_a_line(fp, line, STRLEN))
 115     {
 116         i = line_no;
 117         if (sscanf(line, "%s %d %d %lf", atomnm, &p, &n, &slength) == 4)
 118         {
 119             gnsf->atomnm[i]  = strdup(atomnm);
 120             gnsf->n[i]       = n;
 121             gnsf->p[i]       = p;
 122             gnsf->slength[i] = slength;
 123             line_no++;
 124             gnsf->nratoms = line_no;
 125             if (line_no == nralloc)
 126             {
 127                 nralloc++;
 128                 srenew(gnsf->atomnm, nralloc);
 129                 srenew(gnsf->p, nralloc);
 130                 srenew(gnsf->n, nralloc);
 131                 srenew(gnsf->slength, nralloc);
 132             }
 133         }
 134         else
 135         {
 136             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 137                     datfn, line_no);
 138         }
 139     }
 140     srenew(gnsf->atomnm, gnsf->nratoms);
 141     srenew(gnsf->p, gnsf->nratoms);
 142     srenew(gnsf->n, gnsf->nratoms);
 143     srenew(gnsf->slength, gnsf->nratoms);
 144
 145     fclose(fp);
 146
 147     return (gmx_neutron_atomic_structurefactors_t *) gnsf;
 148 }
 149
 150 gmx_sans_t *gmx_sans_init (t_topology *top, gmx_neutron_atomic_structurefactors_t *gnsf)
 151 {
 152     gmx_sans_t    *gsans = NULL;
 153     int            i, j;
 154     /* Try to assing scattering length from nsfactor.dat */
 155     snew(gsans, 1);
 156     snew(gsans->slength, top->atoms.nr);
 157     /* copy topology data */
 158     gsans->top = top;
 159     for (i = 0; i < top->atoms.nr; i++)
 160     {
 161         for (j = 0; j < gnsf->nratoms; j++)
 162         {
 163             if (top->atoms.atom[i].atomnumber == gnsf->p[j])
 164             {
 165                 /* we need special case for H and D */
 166                 if (top->atoms.atom[i].atomnumber == 1)
 167                 {
 168                     if (top->atoms.atom[i].m == 1.008000)
 169                     {
 170                         gsans->slength[i] = gnsf->slength[0];
 171                     }
 172                     else
 173                     {
 174                         gsans->slength[i] = gnsf->slength[1];
 175                     }
 176                 }
 177                 else
 178                 {
 179                     gsans->slength[i] = gnsf->slength[j];
 180                 }
 181             }
 182         }
 183     }
 184
 185     return (gmx_sans_t *) gsans;
 186 }
 187
 188 gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
 189         gmx_sans_t  *gsans,
 190         rvec        *x,
 191         matrix       box,
 192         atom_id     *index,
 193         int          isize,
 194         double       binwidth,
 195         gmx_bool     bMC,
 196         gmx_bool     bNORM,
 197         real         mcover,
 198         unsigned int seed)
 199 {
 200     gmx_radial_distribution_histogram_t    *pr = NULL;
 201     rvec              dist;
 202     double            rmax;
 203     int               i, j;
 204 #ifdef GMX_OPENMP
 205     double          **tgr;
 206     int               tid;
 207     int               nthreads;
 208     gmx_rng_t        *trng = NULL;
 209 #endif
 210     gmx_int64_t       mc  = 0, max;
 211     gmx_rng_t         rng = NULL;
 212
 213     /* allocate memory for pr */
 214     snew(pr, 1);
 215     /* set some fields */
 216     pr->binwidth = binwidth;
 217
 218     /*
 219      * create max dist rvec
 220      * dist = box[xx] + box[yy] + box[zz]
 221      */
 222     rvec_add(box[XX], box[YY], dist);
 223     rvec_add(box[ZZ], dist, dist);
 224
 225     rmax = norm(dist);
 226
 227     pr->grn = (int)floor(rmax/pr->binwidth)+1;
 228     rmax    = pr->grn*pr->binwidth;
 229
 230     snew(pr->gr, pr->grn);
 231
 232     if (bMC)
 233     {
 234         /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
 235         if (mcover == -1)
 236         {
 237             max = (gmx_int64_t)floor(0.5*0.01*isize*(isize-1));
 238         }
 239         else
 240         {
 241             max = (gmx_int64_t)floor(0.5*mcover*isize*(isize-1));
 242         }
 243         rng = gmx_rng_init(seed);
 244 #ifdef GMX_OPENMP
 245         nthreads = gmx_omp_get_max_threads();
 246         snew(tgr, nthreads);
 247         snew(trng, nthreads);
 248         for (i = 0; i < nthreads; i++)
 249         {
 250             snew(tgr[i], pr->grn);
 251             trng[i] = gmx_rng_init(gmx_rng_uniform_uint32(rng));
 252         }
 253         #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
 254         {
 255             tid = gmx_omp_get_thread_num();
 256             /* now starting parallel threads */
 257             #pragma omp for
 258             for (mc = 0; mc < max; mc++)
 259             {
 260                 i = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
 261                 j = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
 262                 if (i != j)
 263                 {
 264                     tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 265                 }
 266             }
 267         }
 268         /* collecting data from threads */
 269         for (i = 0; i < pr->grn; i++)
 270         {
 271             for (j = 0; j < nthreads; j++)
 272             {
 273                 pr->gr[i] += tgr[j][i];
 274             }
 275         }
 276         /* freeing memory for tgr and destroying trng */
 277         for (i = 0; i < nthreads; i++)
 278         {
 279             sfree(tgr[i]);
 280             gmx_rng_destroy(trng[i]);
 281         }
 282         sfree(tgr);
 283         sfree(trng);
 284 #else
 285         for (mc = 0; mc < max; mc++)
 286         {
 287             i = (int)floor(gmx_rng_uniform_real(rng)*isize);
 288             j = (int)floor(gmx_rng_uniform_real(rng)*isize);
 289             if (i != j)
 290             {
 291                 pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 292             }
 293         }
 294 #endif
 295         gmx_rng_destroy(rng);
 296     }
 297     else
 298     {
 299 #ifdef GMX_OPENMP
 300         nthreads = gmx_omp_get_max_threads();
 301         /* Allocating memory for tgr arrays */
 302         snew(tgr, nthreads);
 303         for (i = 0; i < nthreads; i++)
 304         {
 305             snew(tgr[i], pr->grn);
 306         }
 307         #pragma omp parallel shared(tgr) private(tid,i,j)
 308         {
 309             tid = gmx_omp_get_thread_num();
 310             /* starting parallel threads */
 311             #pragma omp for
 312             for (i = 0; i < isize; i++)
 313             {
 314                 for (j = 0; j < i; j++)
 315                 {
 316                     tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 317                 }
 318             }
 319         }
 320         /* collecating data for pr->gr */
 321         for (i = 0; i < pr->grn; i++)
 322         {
 323             for (j = 0; j < nthreads; j++)
 324             {
 325                 pr->gr[i] += tgr[j][i];
 326             }
 327         }
 328         /* freeing memory for tgr */
 329         for (i = 0; i < nthreads; i++)
 330         {
 331             sfree(tgr[i]);
 332         }
 333         sfree(tgr);
 334 #else
 335         for (i = 0; i < isize; i++)
 336         {
 337             for (j = 0; j < i; j++)
 338             {
 339                 pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 340             }
 341         }
 342 #endif
 343     }
 344
 345     /* normalize if needed */
 346     if (bNORM)
 347     {
 348         normalize_probability(pr->grn, pr->gr);
 349     }
 350
 351     snew(pr->r, pr->grn);
 352     for (i = 0; i < pr->grn; i++)
 353     {
 354         pr->r[i] = (pr->binwidth*i+pr->binwidth*0.5);
 355     }
 356
 357     return (gmx_radial_distribution_histogram_t *) pr;
 358 }
 359
 360 gmx_static_structurefactor_t *convert_histogram_to_intensity_curve (gmx_radial_distribution_histogram_t *pr, double start_q, double end_q, double q_step)
 361 {
 362     gmx_static_structurefactor_t    *sq = NULL;
 363     int         i, j;
 364     /* init data */
 365     snew(sq, 1);
 366     sq->qn = (int)floor((end_q-start_q)/q_step);
 367     snew(sq->q, sq->qn);
 368     snew(sq->s, sq->qn);
 369     for (i = 0; i < sq->qn; i++)
 370     {
 371         sq->q[i] = start_q+i*q_step;
 372     }
 373
 374     if (start_q == 0.0)
 375     {
 376         sq->s[0] = 1.0;
 377         for (i = 1; i < sq->qn; i++)
 378         {
 379             for (j = 0; j < pr->grn; j++)
 380             {
 381                 sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
 382             }
 383             sq->s[i] /= sq->q[i];
 384         }
 385     }
 386     else
 387     {
 388         for (i = 0; i < sq->qn; i++)
 389         {
 390             for (j = 0; j < pr->grn; j++)
 391             {
 392                 sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
 393             }
 394             sq->s[i] /= sq->q[i];
 395         }
 396     }
 397
 398     return (gmx_static_structurefactor_t *) sq;
 399 }