src/gromacs/gmxana/nsfactor.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015, 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 <cmath>
  42 #include <cstring>
  43
  44 #include "gromacs/fileio/strdb.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/random/random.h"
  47 #include "gromacs/topology/topology.h"
  48 #include "gromacs/utility/cstringutil.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/futil.h"
  51 #include "gromacs/utility/gmxomp.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 void check_binwidth(real binwidth)
  55 {
  56     real smallest_bin = 0.1;
  57     if (binwidth < smallest_bin)
  58     {
  59         gmx_fatal(FARGS, "Binwidth shouldnt be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
  60     }
  61 }
  62
  63 void check_mcover(real mcover)
  64 {
  65     if (mcover > 1.0)
  66     {
  67         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  68     }
  69     else if ((mcover < 0)&(mcover != -1))
  70     {
  71         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  72     }
  73     else
  74     {
  75         return;
  76     }
  77 }
  78
  79 void normalize_probability(int n, double *a)
  80 {
  81     int    i;
  82     double norm = 0.0;
  83     for (i = 0; i < n; i++)
  84     {
  85         norm += a[i];
  86     }
  87     for (i = 0; i < n; i++)
  88     {
  89         a[i] /= norm;
  90     }
  91 }
  92
  93 gmx_neutron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn)
  94 {
  95     /* read nsfactor.dat */
  96     FILE    *fp;
  97     char     line[STRLEN];
  98     int      nralloc = 10;
  99     int      n, p;
 100     int      i, line_no;
 101     char     atomnm[8];
 102     double   slength;
 103     gmx_neutron_atomic_structurefactors_t   *gnsf;
 104
 105     fp      = libopen(datfn);
 106     line_no = 0;
 107     /* allocate memory for structure */
 108     snew(gnsf, nralloc);
 109     snew(gnsf->atomnm, nralloc);
 110     snew(gnsf->p, nralloc);
 111     snew(gnsf->n, nralloc);
 112     snew(gnsf->slength, nralloc);
 113
 114     gnsf->nratoms = line_no;
 115
 116     while (get_a_line(fp, line, STRLEN))
 117     {
 118         i = line_no;
 119         if (sscanf(line, "%s %d %d %lf", atomnm, &p, &n, &slength) == 4)
 120         {
 121             gnsf->atomnm[i]  = gmx_strdup(atomnm);
 122             gnsf->n[i]       = n;
 123             gnsf->p[i]       = p;
 124             gnsf->slength[i] = slength;
 125             line_no++;
 126             gnsf->nratoms = line_no;
 127             if (line_no == nralloc)
 128             {
 129                 nralloc++;
 130                 srenew(gnsf->atomnm, nralloc);
 131                 srenew(gnsf->p, nralloc);
 132                 srenew(gnsf->n, nralloc);
 133                 srenew(gnsf->slength, nralloc);
 134             }
 135         }
 136         else
 137         {
 138             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 139                     datfn, line_no);
 140         }
 141     }
 142     srenew(gnsf->atomnm, gnsf->nratoms);
 143     srenew(gnsf->p, gnsf->nratoms);
 144     srenew(gnsf->n, gnsf->nratoms);
 145     srenew(gnsf->slength, gnsf->nratoms);
 146
 147     fclose(fp);
 148
 149     return (gmx_neutron_atomic_structurefactors_t *) gnsf;
 150 }
 151
 152 gmx_sans_t *gmx_sans_init (t_topology *top, gmx_neutron_atomic_structurefactors_t *gnsf)
 153 {
 154     gmx_sans_t    *gsans = NULL;
 155     int            i, j;
 156     /* Try to assing scattering length from nsfactor.dat */
 157     snew(gsans, 1);
 158     snew(gsans->slength, top->atoms.nr);
 159     /* copy topology data */
 160     gsans->top = top;
 161     for (i = 0; i < top->atoms.nr; i++)
 162     {
 163         for (j = 0; j < gnsf->nratoms; j++)
 164         {
 165             if (top->atoms.atom[i].atomnumber == gnsf->p[j])
 166             {
 167                 /* we need special case for H and D */
 168                 if (top->atoms.atom[i].atomnumber == 1)
 169                 {
 170                     if (top->atoms.atom[i].m == 1.008000)
 171                     {
 172                         gsans->slength[i] = gnsf->slength[0];
 173                     }
 174                     else
 175                     {
 176                         gsans->slength[i] = gnsf->slength[1];
 177                     }
 178                 }
 179                 else
 180                 {
 181                     gsans->slength[i] = gnsf->slength[j];
 182                 }
 183             }
 184         }
 185     }
 186
 187     return (gmx_sans_t *) gsans;
 188 }
 189
 190 gmx_radial_distribution_histogram_t *calc_radial_distribution_histogram (
 191         gmx_sans_t  *gsans,
 192         rvec        *x,
 193         matrix       box,
 194         atom_id     *index,
 195         int          isize,
 196         double       binwidth,
 197         gmx_bool     bMC,
 198         gmx_bool     bNORM,
 199         real         mcover,
 200         unsigned int seed)
 201 {
 202     gmx_radial_distribution_histogram_t    *pr = NULL;
 203     rvec              dist;
 204     double            rmax;
 205     int               i, j;
 206 #ifdef GMX_OPENMP
 207     double          **tgr;
 208     int               tid;
 209     int               nthreads;
 210     gmx_rng_t        *trng = NULL;
 211 #endif
 212     gmx_int64_t       mc  = 0, mc_max;
 213     gmx_rng_t         rng = NULL;
 214
 215     /* allocate memory for pr */
 216     snew(pr, 1);
 217     /* set some fields */
 218     pr->binwidth = binwidth;
 219
 220     /*
 221      * create max dist rvec
 222      * dist = box[xx] + box[yy] + box[zz]
 223      */
 224     rvec_add(box[XX], box[YY], dist);
 225     rvec_add(box[ZZ], dist, dist);
 226
 227     rmax = norm(dist);
 228
 229     pr->grn = static_cast<int>(std::floor(rmax/pr->binwidth)+1);
 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             mc_max = static_cast<gmx_int64_t>(std::floor(0.5*0.01*isize*(isize-1)));
 239         }
 240         else
 241         {
 242             mc_max = static_cast<gmx_int64_t>(std::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 < mc_max; mc++)
 260             {
 261                 i = static_cast<int>(std::floor(gmx_rng_uniform_real(trng[tid])*isize));
 262                 j = static_cast<int>(std::floor(gmx_rng_uniform_real(trng[tid])*isize));
 263                 if (i != j)
 264                 {
 265                     tgr[tid][static_cast<int>(std::floor(std::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 < mc_max; mc++)
 287         {
 288             i = static_cast<int>(std::floor(gmx_rng_uniform_real(rng)*isize));
 289             j = static_cast<int>(std::floor(gmx_rng_uniform_real(rng)*isize));
 290             if (i != j)
 291             {
 292                 pr->gr[static_cast<int>(std::floor(std::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][static_cast<int>(std::floor(std::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[static_cast<int>(std::floor(std::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 = static_cast<int>(std::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 }