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,2016,2017,2018, 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/math/vec.h"
  45 #include "gromacs/random/threefry.h"
  46 #include "gromacs/random/uniformintdistribution.h"
  47 #include "gromacs/topology/topology.h"
  48 #include "gromacs/utility/cstringutil.h"
  49 #include "gromacs/utility/exceptions.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/futil.h"
  52 #include "gromacs/utility/gmxomp.h"
  53 #include "gromacs/utility/smalloc.h"
  54 #include "gromacs/utility/strdb.h"
  55
  56 void check_binwidth(real binwidth)
  57 {
  58     real smallest_bin = 0.1;
  59     if (binwidth < smallest_bin)
  60     {
  61         gmx_fatal(FARGS, "Binwidth shouldn't be smaller then smallest bond length (H-H bond ~0.1nm) in a box");
  62     }
  63 }
  64
  65 void check_mcover(real mcover)
  66 {
  67     if (mcover > 1.0)
  68     {
  69         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  70     }
  71     else if ((mcover < 0) && (mcover != -1))
  72     {
  73         gmx_fatal(FARGS, "mcover should be -1 or (0,1]");
  74     }
  75     else
  76     {
  77         return;
  78     }
  79 }
  80
  81 void normalize_probability(int n, double *a)
  82 {
  83     int    i;
  84     double norm = 0.0;
  85     for (i = 0; i < n; i++)
  86     {
  87         norm += a[i];
  88     }
  89     for (i = 0; i < n; i++)
  90     {
  91         a[i] /= norm;
  92     }
  93 }
  94
  95 gmx_neutron_atomic_structurefactors_t *gmx_neutronstructurefactors_init(const char *datfn)
  96 {
  97     /* read nsfactor.dat */
  98     char     line[STRLEN];
  99     int      nralloc = 10;
 100     int      n, p;
 101     int      i, line_no;
 102     char     atomnm[8];
 103     double   slength;
 104     gmx_neutron_atomic_structurefactors_t   *gnsf;
 105
 106     gmx::FilePtr fp = gmx::openLibraryFile(datfn);
 107     line_no = 0;
 108     /* allocate memory for structure */
 109     snew(gnsf, nralloc);
 110     snew(gnsf->atomnm, nralloc);
 111     snew(gnsf->p, nralloc);
 112     snew(gnsf->n, nralloc);
 113     snew(gnsf->slength, nralloc);
 114
 115     gnsf->nratoms = line_no;
 116
 117     while (get_a_line(fp.get(), line, STRLEN))
 118     {
 119         i = line_no;
 120         if (sscanf(line, "%s %d %d %lf", atomnm, &p, &n, &slength) == 4)
 121         {
 122             gnsf->atomnm[i]  = gmx_strdup(atomnm);
 123             gnsf->n[i]       = n;
 124             gnsf->p[i]       = p;
 125             gnsf->slength[i] = slength;
 126             line_no++;
 127             gnsf->nratoms = line_no;
 128             if (line_no == nralloc)
 129             {
 130                 nralloc++;
 131                 srenew(gnsf->atomnm, nralloc);
 132                 srenew(gnsf->p, nralloc);
 133                 srenew(gnsf->n, nralloc);
 134                 srenew(gnsf->slength, nralloc);
 135             }
 136         }
 137         else
 138         {
 139             fprintf(stderr, "WARNING: Error in file %s at line %d ignored\n",
 140                     datfn, line_no);
 141         }
 142     }
 143     srenew(gnsf->atomnm, gnsf->nratoms);
 144     srenew(gnsf->p, gnsf->nratoms);
 145     srenew(gnsf->n, gnsf->nratoms);
 146     srenew(gnsf->slength, gnsf->nratoms);
 147
 148     return gnsf;
 149 }
 150
 151 gmx_sans_t *gmx_sans_init (const t_topology *top, gmx_neutron_atomic_structurefactors_t *gnsf)
 152 {
 153     gmx_sans_t    *gsans = nullptr;
 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 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         const int         *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 = nullptr;
 202     rvec                                    dist;
 203     double                                  rmax;
 204     int                                     i, j;
 205 #if GMX_OPENMP
 206     double                                **tgr;
 207     int                                     tid;
 208     int                                     nthreads;
 209     gmx::DefaultRandomEngine               *trng = nullptr;
 210 #endif
 211     int64_t                                 mc  = 0, mc_max;
 212     gmx::DefaultRandomEngine                rng(seed);
 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 = static_cast<int>(std::floor(rmax/pr->binwidth)+1);
 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             mc_max = static_cast<int64_t>(std::floor(0.5*0.01*isize*(isize-1)));
 238         }
 239         else
 240         {
 241             mc_max = static_cast<int64_t>(std::floor(0.5*mcover*isize*(isize-1)));
 242         }
 243 #if GMX_OPENMP
 244         nthreads = gmx_omp_get_max_threads();
 245         snew(tgr, nthreads);
 246         trng = new gmx::DefaultRandomEngine[nthreads];
 247         for (i = 0; i < nthreads; i++)
 248         {
 249             snew(tgr[i], pr->grn);
 250             trng[i].seed(rng());
 251         }
 252         #pragma omp parallel shared(tgr,trng,mc) private(tid,i,j)
 253         {
 254             gmx::UniformIntDistribution<int> tdist(0, isize-1);
 255             tid = gmx_omp_get_thread_num();
 256             /* now starting parallel threads */
 257             #pragma omp for
 258             for (mc = 0; mc < mc_max; mc++)
 259             {
 260                 try
 261                 {
 262                     i = tdist(trng[tid]); // [0,isize-1]
 263                     j = tdist(trng[tid]); // [0,isize-1]
 264                     if (i != j)
 265                     {
 266                         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]];
 267                     }
 268                 }
 269                 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 270             }
 271         }
 272         /* collecting data from threads */
 273         for (i = 0; i < pr->grn; i++)
 274         {
 275             for (j = 0; j < nthreads; j++)
 276             {
 277                 pr->gr[i] += tgr[j][i];
 278             }
 279         }
 280         /* freeing memory for tgr and destroying trng */
 281         for (i = 0; i < nthreads; i++)
 282         {
 283             sfree(tgr[i]);
 284         }
 285         sfree(tgr);
 286         delete[] trng;
 287 #else
 288         gmx::UniformIntDistribution<int> dist(0, isize-1);
 289         for (mc = 0; mc < mc_max; mc++)
 290         {
 291             i = dist(rng); // [0,isize-1]
 292             j = dist(rng); // [0,isize-1]
 293             if (i != j)
 294             {
 295                 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]];
 296             }
 297         }
 298 #endif
 299     }
 300     else
 301     {
 302 #if GMX_OPENMP
 303         nthreads = gmx_omp_get_max_threads();
 304         /* Allocating memory for tgr arrays */
 305         snew(tgr, nthreads);
 306         for (i = 0; i < nthreads; i++)
 307         {
 308             snew(tgr[i], pr->grn);
 309         }
 310         #pragma omp parallel shared(tgr) private(tid,i,j)
 311         {
 312             tid = gmx_omp_get_thread_num();
 313             /* starting parallel threads */
 314             #pragma omp for
 315             for (i = 0; i < isize; i++)
 316             {
 317                 try
 318                 {
 319                     for (j = 0; j < i; j++)
 320                     {
 321                         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]];
 322                     }
 323                 }
 324                 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
 325             }
 326         }
 327         /* collecating data for pr->gr */
 328         for (i = 0; i < pr->grn; i++)
 329         {
 330             for (j = 0; j < nthreads; j++)
 331             {
 332                 pr->gr[i] += tgr[j][i];
 333             }
 334         }
 335         /* freeing memory for tgr */
 336         for (i = 0; i < nthreads; i++)
 337         {
 338             sfree(tgr[i]);
 339         }
 340         sfree(tgr);
 341 #else
 342         for (i = 0; i < isize; i++)
 343         {
 344             for (j = 0; j < i; j++)
 345             {
 346                 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]];
 347             }
 348         }
 349 #endif
 350     }
 351
 352     /* normalize if needed */
 353     if (bNORM)
 354     {
 355         normalize_probability(pr->grn, pr->gr);
 356     }
 357
 358     snew(pr->r, pr->grn);
 359     for (i = 0; i < pr->grn; i++)
 360     {
 361         pr->r[i] = (pr->binwidth*i+pr->binwidth*0.5);
 362     }
 363
 364     return pr;
 365 }
 366
 367 gmx_static_structurefactor_t *convert_histogram_to_intensity_curve (gmx_radial_distribution_histogram_t *pr, double start_q, double end_q, double q_step)
 368 {
 369     gmx_static_structurefactor_t    *sq = nullptr;
 370     int         i, j;
 371     /* init data */
 372     snew(sq, 1);
 373     sq->qn = static_cast<int>(std::floor((end_q-start_q)/q_step));
 374     snew(sq->q, sq->qn);
 375     snew(sq->s, sq->qn);
 376     for (i = 0; i < sq->qn; i++)
 377     {
 378         sq->q[i] = start_q+i*q_step;
 379     }
 380
 381     if (start_q == 0.0)
 382     {
 383         sq->s[0] = 1.0;
 384         for (i = 1; i < sq->qn; i++)
 385         {
 386             for (j = 0; j < pr->grn; j++)
 387             {
 388                 sq->s[i] += (pr->gr[j]/pr->r[j])*std::sin(sq->q[i]*pr->r[j]);
 389             }
 390             sq->s[i] /= sq->q[i];
 391         }
 392     }
 393     else
 394     {
 395         for (i = 0; i < sq->qn; i++)
 396         {
 397             for (j = 0; j < pr->grn; j++)
 398             {
 399                 sq->s[i] += (pr->gr[j]/pr->r[j])*std::sin(sq->q[i]*pr->r[j]);
 400             }
 401             sq->s[i] /= sq->q[i];
 402         }
 403     }
 404
 405     return sq;
 406 }