src/tools/nsfactor.c

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