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 #include "gmxpre.h"
  36
  37 #include "nsfactor.h"
  38
  39 #include "config.h"
  40
  41 #include <string.h>
  42
  43 #include "gromacs/fileio/strdb.h"
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/random/random.h"
  46 #include "gromacs/topology/topology.h"
  47 #include "gromacs/utility/cstringutil.h"
  48 #include "gromacs/utility/fatalerror.h"
  49 #include "gromacs/utility/futil.h"
  50 #include "gromacs/utility/gmxomp.h"
  51 #include "gromacs/utility/smalloc.h"
  52
  53 void check_binwidth(real binwidth)
  54 {
  55     real smallest_bin = 0.1;
  56     if (binwidth < smallest_bin)
  57     {
  58         gmx_fatal(FARGS, "Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
  59     }
  60 }
  61
  62 void check_mcover(real mcover)
  63 {
  64     if (mcover > 1.0)
  65     {
  66         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  67     }
  68     else if ((mcover < 0)&(mcover != -1))
  69     {
  70         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  71     }
  72     else
  73     {
  74         return;
  75     }
  76 }
  77
  78 void normalize_probability(int n, double *a)
  79 {
  80     int    i;
  81     double norm = 0.0;
  82     for (i = 0; i < n; i++)
  83     {
  84         norm += a[i];
  85     }
  86     for (i = 0; i < n; i++)
  87     {
  88         a[i] /= norm;
  89     }
  90 }
  91
  92 gmx_neutron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn)
  93 {
  94     /* read nsfactor.dat */
  95     FILE    *fp;
  96     char     line[STRLEN];
  97     int      nralloc = 10;
  98     int      n, p;
  99     int      i, line_no;
 100     char     atomnm[8];
 101     double   slength;
 102     gmx_neutron_atomic_structurefactors_t   *gnsf;
 103
 104     fp      = libopen(datfn);
 105     line_no = 0;
 106     /* allocate memory for structure */
 107     snew(gnsf, nralloc);
 108     snew(gnsf->atomnm, nralloc);
 109     snew(gnsf->p, nralloc);
 110     snew(gnsf->n, nralloc);
 111     snew(gnsf->slength, nralloc);
 112
 113     gnsf->nratoms = line_no;
 114
 115     while (get_a_line(fp, line, STRLEN))
 116     {
 117         i = line_no;
 118         if (sscanf(line, "%s %d %d %lf", atomnm, &p, &n, &slength) == 4)
 119         {
 120             gnsf->atomnm[i]  = gmx_strdup(atomnm);
 121             gnsf->n[i]       = n;
 122             gnsf->p[i]       = p;
 123             gnsf->slength[i] = slength;
 124             line_no++;
 125             gnsf->nratoms = line_no;
 126             if (line_no == nralloc)
 127             {
 128                 nralloc++;
 129                 srenew(gnsf->atomnm, nralloc);
 130                 srenew(gnsf->p, nralloc);
 131                 srenew(gnsf->n, nralloc);
 132                 srenew(gnsf->slength, nralloc);
 133             }
 134         }
 135         else
 136         {
 137             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 138                     datfn, line_no);
 139         }
 140     }
 141     srenew(gnsf->atomnm, gnsf->nratoms);
 142     srenew(gnsf->p, gnsf->nratoms);
 143     srenew(gnsf->n, gnsf->nratoms);
 144     srenew(gnsf->slength, gnsf->nratoms);
 145
 146     fclose(fp);
 147
 148     return (gmx_neutron_atomic_structurefactors_t *) gnsf;
 149 }
 150
 151 gmx_sans_t *gmx_sans_init (t_topology *top, gmx_neutron_atomic_structurefactors_t *gnsf)
 152 {
 153     gmx_sans_t    *gsans = NULL;
 154     int            i, j;
 155     /* Try to assing scattering length from nsfactor.dat */
 156     snew(gsans, 1);
 157     snew(gsans->slength, top->atoms.nr);
 158     /* copy topology data */
 159     gsans->top = top;
 160     for (i = 0; i < top->atoms.nr; i++)
 161     {
 162         for (j = 0; j < gnsf->nratoms; j++)
 163         {
 164             if (top->atoms.atom[i].atomnumber == gnsf->p[j])
 165             {
 166                 /* we need special case for H and D */
 167                 if (top->atoms.atom[i].atomnumber == 1)
 168                 {
 169                     if (top->atoms.atom[i].m == 1.008000)
 170                     {
 171                         gsans->slength[i] = gnsf->slength[0];
 172                     }
 173                     else
 174                     {
 175                         gsans->slength[i] = gnsf->slength[1];
 176                     }
 177                 }
 178                 else
 179                 {
 180                     gsans->slength[i] = gnsf->slength[j];
 181                 }
 182             }
 183         }
 184     }
 185
 186     return (gmx_sans_t *) gsans;
 187 }
 188
 189 gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
 190         gmx_sans_t  *gsans,
 191         rvec        *x,
 192         matrix       box,
 193         atom_id     *index,
 194         int          isize,
 195         double       binwidth,
 196         gmx_bool     bMC,
 197         gmx_bool     bNORM,
 198         real         mcover,
 199         unsigned int seed)
 200 {
 201     gmx_radial_distribution_histogram_t    *pr = NULL;
 202     rvec              dist;
 203     double            rmax;
 204     int               i, j;
 205 #ifdef GMX_OPENMP
 206     double          **tgr;
 207     int               tid;
 208     int               nthreads;
 209     gmx_rng_t        *trng = NULL;
 210 #endif
 211     gmx_int64_t       mc  = 0, max;
 212     gmx_rng_t         rng = NULL;
 213
 214     /* allocate memory for pr */
 215     snew(pr, 1);
 216     /* set some fields */
 217     pr->binwidth = binwidth;
 218
 219     /*
 220      * create max dist rvec
 221      * dist = box[xx] + box[yy] + box[zz]
 222      */
 223     rvec_add(box[XX], box[YY], dist);
 224     rvec_add(box[ZZ], dist, dist);
 225
 226     rmax = norm(dist);
 227
 228     pr->grn = (int)floor(rmax/pr->binwidth)+1;
 229     rmax    = pr->grn*pr->binwidth;
 230
 231     snew(pr->gr, pr->grn);
 232
 233     if (bMC)
 234     {
 235         /* Special case for setting automaticaly number of mc iterations to 1% of total number of direct iterations */
 236         if (mcover == -1)
 237         {
 238             max = (gmx_int64_t)floor(0.5*0.01*isize*(isize-1));
 239         }
 240         else
 241         {
 242             max = (gmx_int64_t)floor(0.5*mcover*isize*(isize-1));
 243         }
 244         rng = gmx_rng_init(seed);
 245 #ifdef GMX_OPENMP
 246         nthreads = gmx_omp_get_max_threads();
 247         snew(tgr, nthreads);
 248         snew(trng, nthreads);
 249         for (i = 0; i < nthreads; i++)
 250         {
 251             snew(tgr[i], pr->grn);
 252             trng[i] = gmx_rng_init(gmx_rng_uniform_uint32(rng));
 253         }
 254         #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
 255         {
 256             tid = gmx_omp_get_thread_num();
 257             /* now starting parallel threads */
 258             #pragma omp for
 259             for (mc = 0; mc < max; mc++)
 260             {
 261                 i = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
 262                 j = (int)floor(gmx_rng_uniform_real(trng[tid])*isize);
 263                 if (i != j)
 264                 {
 265                     tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 266                 }
 267             }
 268         }
 269         /* collecting data from threads */
 270         for (i = 0; i < pr->grn; i++)
 271         {
 272             for (j = 0; j < nthreads; j++)
 273             {
 274                 pr->gr[i] += tgr[j][i];
 275             }
 276         }
 277         /* freeing memory for tgr and destroying trng */
 278         for (i = 0; i < nthreads; i++)
 279         {
 280             sfree(tgr[i]);
 281             gmx_rng_destroy(trng[i]);
 282         }
 283         sfree(tgr);
 284         sfree(trng);
 285 #else
 286         for (mc = 0; mc < max; mc++)
 287         {
 288             i = (int)floor(gmx_rng_uniform_real(rng)*isize);
 289             j = (int)floor(gmx_rng_uniform_real(rng)*isize);
 290             if (i != j)
 291             {
 292                 pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 293             }
 294         }
 295 #endif
 296         gmx_rng_destroy(rng);
 297     }
 298     else
 299     {
 300 #ifdef GMX_OPENMP
 301         nthreads = gmx_omp_get_max_threads();
 302         /* Allocating memory for tgr arrays */
 303         snew(tgr, nthreads);
 304         for (i = 0; i < nthreads; i++)
 305         {
 306             snew(tgr[i], pr->grn);
 307         }
 308         #pragma omp parallel shared(tgr) private(tid,i,j)
 309         {
 310             tid = gmx_omp_get_thread_num();
 311             /* starting parallel threads */
 312             #pragma omp for
 313             for (i = 0; i < isize; i++)
 314             {
 315                 for (j = 0; j < i; j++)
 316                 {
 317                     tgr[tid][(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 318                 }
 319             }
 320         }
 321         /* collecating data for pr->gr */
 322         for (i = 0; i < pr->grn; i++)
 323         {
 324             for (j = 0; j < nthreads; j++)
 325             {
 326                 pr->gr[i] += tgr[j][i];
 327             }
 328         }
 329         /* freeing memory for tgr */
 330         for (i = 0; i < nthreads; i++)
 331         {
 332             sfree(tgr[i]);
 333         }
 334         sfree(tgr);
 335 #else
 336         for (i = 0; i < isize; i++)
 337         {
 338             for (j = 0; j < i; j++)
 339             {
 340                 pr->gr[(int)floor(sqrt(distance2(x[index[i]], x[index[j]]))/binwidth)] += gsans->slength[index[i]]*gsans->slength[index[j]];
 341             }
 342         }
 343 #endif
 344     }
 345
 346     /* normalize if needed */
 347     if (bNORM)
 348     {
 349         normalize_probability(pr->grn, pr->gr);
 350     }
 351
 352     snew(pr->r, pr->grn);
 353     for (i = 0; i < pr->grn; i++)
 354     {
 355         pr->r[i] = (pr->binwidth*i+pr->binwidth*0.5);
 356     }
 357
 358     return (gmx_radial_distribution_histogram_t *) pr;
 359 }
 360
 361 gmx_static_structurefactor_t *convert_histogram_to_intensity_curve (gmx_radial_distribution_histogram_t *pr, double start_q, double end_q, double q_step)
 362 {
 363     gmx_static_structurefactor_t    *sq = NULL;
 364     int         i, j;
 365     /* init data */
 366     snew(sq, 1);
 367     sq->qn = (int)floor((end_q-start_q)/q_step);
 368     snew(sq->q, sq->qn);
 369     snew(sq->s, sq->qn);
 370     for (i = 0; i < sq->qn; i++)
 371     {
 372         sq->q[i] = start_q+i*q_step;
 373     }
 374
 375     if (start_q == 0.0)
 376     {
 377         sq->s[0] = 1.0;
 378         for (i = 1; i < sq->qn; i++)
 379         {
 380             for (j = 0; j < pr->grn; j++)
 381             {
 382                 sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
 383             }
 384             sq->s[i] /= sq->q[i];
 385         }
 386     }
 387     else
 388     {
 389         for (i = 0; i < sq->qn; i++)
 390         {
 391             for (j = 0; j < pr->grn; j++)
 392             {
 393                 sq->s[i] += (pr->gr[j]/pr->r[j])*sin(sq->q[i]*pr->r[j]);
 394             }
 395             sq->s[i] /= sq->q[i];
 396         }
 397     }
 398
 399     return (gmx_static_structurefactor_t *) sq;
 400 }