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