src/tools/gmx_sans.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <ctype.h>
  44 #include "smalloc.h"
  45 #include "sysstuff.h"
  46 #include "typedefs.h"
  47 #include "macros.h"
  48 #include "vec.h"
  49 #include "pbc.h"
  50 #include "xvgr.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "statutil.h"
  54 #include "tpxio.h"
  55 #include "index.h"
  56 #include "gstat.h"
  57 #include "matio.h"
  58 #include "gmx_ana.h"
  59 #include "nsfactor.h"
  60 #include "gmx_omp.h"
  61
  62 int gmx_sans(int argc,char *argv[])
  63 {
  64     const char *desc[] = {
  65         "This is simple tool to compute SANS spectra using Debye formula",
  66         "It currently uses topology file (since it need to assigne element for each atom)",
  67         "[PAR]",
  68         "[TT]-pr[tt] Computes normalized g(r) function",
  69         "[PAR]",
  70         "[TT]-sq[tt] Computes SANS intensity curve for needed q diapason",
  71         "[PAR]",
  72         "[TT]-startq[tt] Starting q value in nm",
  73         "[PAR]",
  74         "[TT]-endq[tt] Ending q value in nm",
  75         "[PAR]",
  76         "[TT]-qstep[tt] Stepping in q space",
  77         "[PAR]",
  78         "Note: When using Debye direct method computational cost increases as",
  79         "1/2 * N * (N - 1) where N is atom number in group of interest"
  80     };
  81     static gmx_bool bPBC=TRUE;
  82     static real binwidth=0.2,grid=0.05; /* bins shouldnt be smaller then bond (~0.1nm) length */
  83     static real start_q=0.0, end_q=2.0, q_step=0.01;
  84     static real mcover=-1;
  85     static unsigned int  seed=0;
  86     static int           nthreads=-1;
  87
  88     static const char *emode[]= { NULL, "direct", "mc", NULL };
  89     static const char *emethod[]={ NULL, "debye", "fft", NULL };
  90
  91     gmx_nentron_atomic_structurefactors_t    *gnsf;
  92     gmx_sans_t              *gsans;
  93
  94 #define NPA asize(pa)
  95
  96     t_pargs pa[] = {
  97         { "-bin", FALSE, etREAL, {&binwidth},
  98           "[HIDDEN]Binwidth (nm)" },
  99         { "-mode", FALSE, etENUM, {emode},
 100           "Mode for sans spectra calculation" },
 101         { "-mcover", FALSE, etREAL, {&mcover},
 102           "Monte-Carlo coverage should be -1(default) or (0,1]"},
 103         { "-method", FALSE, etENUM, {emethod},
 104           "[HIDDEN]Method for sans spectra calculation" },
 105         { "-pbc", FALSE, etBOOL, {&bPBC},
 106           "Use periodic boundary conditions for computing distances" },
 107         { "-grid", FALSE, etREAL, {&grid},
 108           "[HIDDEN]Grid spacing (in nm) for FFTs" },
 109         {"-startq", FALSE, etREAL, {&start_q},
 110           "Starting q (1/nm) "},
 111         {"-endq", FALSE, etREAL, {&end_q},
 112           "Ending q (1/nm)"},
 113         { "-qstep", FALSE, etREAL, {&q_step},
 114           "Stepping in q (1/nm)"},
 115         { "-seed",     FALSE, etINT,  {&seed},
 116           "Random seed for Monte-Carlo"},
 117 #ifdef GMX_OPENMP
 118         { "-nt",  FALSE, etINT, {&nthreads},
 119           "Number of threads to start"},
 120 #endif
 121     };
 122   FILE      *fp;
 123   const char *fnTPX,*fnNDX,*fnDAT=NULL;
 124   t_trxstatus *status;
 125   t_topology *top=NULL;
 126   t_atom    *atom=NULL;
 127   gmx_rmpbc_t  gpbc=NULL;
 128   gmx_bool  bTPX;
 129   gmx_bool  bFFT=FALSE, bDEBYE=FALSE;
 130   gmx_bool  bMC=FALSE;
 131   int        ePBC=-1;
 132   matrix     box;
 133   char       title[STRLEN];
 134   rvec       *x;
 135   int       natoms;
 136   real       t;
 137   char       **grpname=NULL;
 138   atom_id    *index=NULL;
 139   int        isize;
 140   int         i,j;
 141   gmx_radial_distribution_histogram_t  *pr=NULL;
 142   gmx_static_structurefator_t  *sq=NULL;
 143   output_env_t oenv;
 144
 145 #define NFILE asize(fnm)
 146
 147   t_filenm   fnm[] = {
 148       { efTPX,  "-s",         NULL,   ffREAD },
 149       { efNDX,  NULL,         NULL,   ffOPTRD },
 150       { efDAT,  "-d",   "nsfactor",   ffOPTRD },
 151       { efXVG, "-sq",         "sq",   ffWRITE },
 152       { efXVG, "-pr",         "pr",   ffWRITE }
 153   };
 154
 155   nthreads = gmx_omp_get_max_threads();
 156
 157   CopyRight(stderr,argv[0]);
 158   parse_common_args(&argc,argv,PCA_BE_NICE,
 159                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 160
 161   /* check that binwidth not smaller than smallers distance */
 162   check_binwidth(binwidth);
 163   check_mcover(mcover);
 164
 165   /* setting number of omp threads globaly */
 166   gmx_omp_set_num_threads(nthreads);
 167
 168   /* Now try to parse opts for modes */
 169   switch(emethod[0][0]) {
 170   case 'd':
 171       bDEBYE=TRUE;
 172       switch(emode[0][0]) {
 173       case 'd':
 174           bMC=FALSE;
 175           break;
 176       case 'm':
 177           bMC=TRUE;
 178           break;
 179       default:
 180           break;
 181       }
 182       break;
 183   case 'f':
 184       bFFT=TRUE;
 185       break;
 186   default:
 187       break;
 188   }
 189
 190   if (!bDEBYE && !bFFT)
 191       gmx_fatal(FARGS,"Unknown method. Set pr or fft!\n");
 192   /* Try to read files */
 193   fnDAT = ftp2fn(efDAT,NFILE,fnm);
 194   fnTPX = ftp2fn(efTPX,NFILE,fnm);
 195
 196   gnsf = gmx_neutronstructurefactors_init(fnDAT);
 197   fprintf(stderr,"Read %d atom names from %s with neutron scattering parameters\n\n",gnsf->nratoms,fnDAT);
 198
 199   snew(top,1);
 200   snew(grpname,1);
 201   snew(index,1);
 202
 203   bTPX=read_tps_conf(fnTPX,title,top,&ePBC,&x,NULL,box,TRUE);
 204
 205   printf("\nPlease select group for SANS spectra calculation:\n");
 206   get_index(&(top->atoms),ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,grpname);
 207
 208   gsans = gmx_sans_init(top,gnsf);
 209
 210   /* Prepare reference frame */
 211   if (bPBC) {
 212       gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
 213       gmx_rmpbc(gpbc,top->atoms.nr,box,x);
 214   }
 215
 216   natoms=top->atoms.nr;
 217
 218   if (bDEBYE) {
 219       if (bMC) {
 220           fprintf(stderr,"Using Monte Carlo Debye method to calculate spectrum\n");
 221       } else {
 222           fprintf(stderr,"Using direct Debye method to calculate spectrum\n");
 223       }
 224   } else if (bFFT) {
 225       gmx_fatal(FARGS,"Not implented!");
 226   } else {
 227       gmx_fatal(FARGS,"Whats this!");
 228   }
 229
 230   /*  realy calc p(r) */
 231   pr = calc_radial_distribution_histogram(gsans,x,box,index,isize,binwidth,bMC,mcover,seed);
 232
 233   /* prepare pr.xvg */
 234   fp = xvgropen(opt2fn_null("-pr",NFILE,fnm),"G(r)","Distance (nm)","Probability",oenv);
 235   for(i=0;i<pr->grn;i++)
 236       fprintf(fp,"%10.6lf%10.6lf\n",pr->r[i],pr->gr[i]);
 237   xvgrclose(fp);
 238
 239   /* prepare sq.xvg */
 240   sq = convert_histogram_to_intensity_curve(pr,start_q,end_q,q_step);
 241   fp = xvgropen(opt2fn_null("-sq",NFILE,fnm),"I(q)","q (nm^-1)","s(q)/s(0)",oenv);
 242   for(i=0;i<sq->qn;i++) {
 243       fprintf(fp,"%10.6lf%10.6lf\n",sq->q[i],sq->s[i]);
 244   }
 245   xvgrclose(fp);
 246
 247   sfree(pr);
 248
 249   please_cite(stdout,"Garmay2012");
 250   thanx(stderr);
 251
 252   return 0;
 253 }