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         "Parameters:[PAR]"
  69         "[TT]-pr[tt] Computes normalized g(r) function averaged over trajectory[PAR]",
  70         "[TT]-prframe[tt] Computes normalized g(r) function for each frame[PAR]",
  71         "[TT]-sq[tt] Computes SANS intensity curve averaged over trajectory[PAR]",
  72         "[TT]-sqframe[tt] Computes SANS intensity curve for each frame[PAR]",
  73         "[TT]-startq[tt] Starting q value in nm[PAR]",
  74         "[TT]-endq[tt] Ending q value in nm[PAR]",
  75         "[TT]-qstep[tt] Stepping in q space[PAR]",
  76         "Note: When using Debye direct method computational cost increases as",
  77         "1/2 * N * (N - 1) where N is atom number in group of interest",
  78         "[PAR]",
  79         "WARNING: If sq or pr specified this tool can produce large number of files! Up to two times larger than number of frames!"
  80     };
  81     static gmx_bool bPBC=TRUE;
  82     static gmx_bool bNORM=FALSE;
  83     static real binwidth=0.2,grid=0.05; /* bins shouldnt be smaller then smallest bond (~0.1nm) length */
  84     static real start_q=0.0, end_q=2.0, q_step=0.01;
  85     static real mcover=-1;
  86     static unsigned int  seed=0;
  87     static int           nthreads=-1;
  88
  89     static const char *emode[]= { NULL, "direct", "mc", NULL };
  90     static const char *emethod[]={ NULL, "debye", "fft", NULL };
  91
  92     gmx_neutron_atomic_structurefactors_t    *gnsf;
  93     gmx_sans_t              *gsans;
  94
  95 #define NPA asize(pa)
  96
  97     t_pargs pa[] = {
  98         { "-bin", FALSE, etREAL, {&binwidth},
  99           "[HIDDEN]Binwidth (nm)" },
 100         { "-mode", FALSE, etENUM, {emode},
 101           "Mode for sans spectra calculation" },
 102         { "-mcover", FALSE, etREAL, {&mcover},
 103           "Monte-Carlo coverage should be -1(default) or (0,1]"},
 104         { "-method", FALSE, etENUM, {emethod},
 105           "[HIDDEN]Method for sans spectra calculation" },
 106         { "-pbc", FALSE, etBOOL, {&bPBC},
 107           "Use periodic boundary conditions for computing distances" },
 108         { "-grid", FALSE, etREAL, {&grid},
 109           "[HIDDEN]Grid spacing (in nm) for FFTs" },
 110         {"-startq", FALSE, etREAL, {&start_q},
 111           "Starting q (1/nm) "},
 112         {"-endq", FALSE, etREAL, {&end_q},
 113           "Ending q (1/nm)"},
 114         { "-qstep", FALSE, etREAL, {&q_step},
 115           "Stepping in q (1/nm)"},
 116         { "-seed",     FALSE, etINT,  {&seed},
 117           "Random seed for Monte-Carlo"},
 118 #ifdef GMX_OPENMP
 119         { "-nt",  FALSE, etINT, {&nthreads},
 120           "Number of threads to start"},
 121 #endif
 122     };
 123   FILE      *fp;
 124   const char *fnTPX,*fnNDX,*fnTRX,*fnDAT=NULL;
 125   t_trxstatus *status;
 126   t_topology *top=NULL;
 127   t_atom    *atom=NULL;
 128   gmx_rmpbc_t  gpbc=NULL;
 129   gmx_bool  bTPX;
 130   gmx_bool  bFFT=FALSE, bDEBYE=FALSE;
 131   gmx_bool  bMC=FALSE;
 132   int        ePBC=-1;
 133   matrix     box;
 134   char       title[STRLEN];
 135   rvec       *x;
 136   int       natoms;
 137   real       t;
 138   char       **grpname=NULL;
 139   atom_id    *index=NULL;
 140   int        isize;
 141   int         i,j;
 142   char       *hdr=NULL;
 143   char       *suffix=NULL;
 144   t_filenm   *fnmdup=NULL;
 145   gmx_radial_distribution_histogram_t  *prframecurrent=NULL, *pr=NULL;
 146   gmx_static_structurefactor_t  *sqframecurrent=NULL, *sq=NULL;
 147   output_env_t oenv;
 148
 149 #define NFILE asize(fnm)
 150
 151   t_filenm   fnm[] = {
 152       { efTPX,  "-s",       NULL,       ffREAD },
 153       { efTRX,  "-f",       NULL,       ffREAD },
 154       { efNDX,  NULL,       NULL,       ffOPTRD },
 155       { efDAT,  "-d",       "nsfactor", ffOPTRD },
 156       { efXVG,  "-pr",      "pr",       ffWRITE },
 157       { efXVG,  "-sq",       "sq",      ffWRITE },
 158       { efXVG,  "-prframe", "prframe",  ffOPTWR },
 159       { efXVG,  "-sqframe", "sqframe",  ffOPTWR }
 160   };
 161
 162   nthreads = gmx_omp_get_max_threads();
 163
 164   CopyRight(stderr,argv[0]);
 165   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 166                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 167
 168   /* check that binwidth not smaller than smallers distance */
 169   check_binwidth(binwidth);
 170   check_mcover(mcover);
 171
 172   /* setting number of omp threads globaly */
 173   gmx_omp_set_num_threads(nthreads);
 174
 175   /* Now try to parse opts for modes */
 176   switch(emethod[0][0]) {
 177   case 'd':
 178       bDEBYE=TRUE;
 179       switch(emode[0][0]) {
 180       case 'd':
 181           bMC=FALSE;
 182           break;
 183       case 'm':
 184           bMC=TRUE;
 185           break;
 186       default:
 187           break;
 188       }
 189       break;
 190   case 'f':
 191       bFFT=TRUE;
 192       break;
 193   default:
 194       break;
 195   }
 196
 197   if (bDEBYE) {
 198       if (bMC) {
 199           fprintf(stderr,"Using Monte Carlo Debye method to calculate spectrum\n");
 200       } else {
 201           fprintf(stderr,"Using direct Debye method to calculate spectrum\n");
 202       }
 203   } else if (bFFT) {
 204       gmx_fatal(FARGS,"FFT method not implemented!");
 205   } else {
 206       gmx_fatal(FARGS,"Unknown combination for mode and method!");
 207   }
 208
 209   /* Try to read files */
 210   fnDAT = ftp2fn(efDAT,NFILE,fnm);
 211   fnTPX = ftp2fn(efTPX,NFILE,fnm);
 212   fnTRX = ftp2fn(efTRX,NFILE,fnm);
 213
 214   gnsf = gmx_neutronstructurefactors_init(fnDAT);
 215   fprintf(stderr,"Read %d atom names from %s with neutron scattering parameters\n\n",gnsf->nratoms,fnDAT);
 216
 217   snew(top,1);
 218   snew(grpname,1);
 219   snew(index,1);
 220
 221   bTPX=read_tps_conf(fnTPX,title,top,&ePBC,&x,NULL,box,TRUE);
 222
 223   printf("\nPlease select group for SANS spectra calculation:\n");
 224   get_index(&(top->atoms),ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,grpname);
 225
 226   gsans = gmx_sans_init(top,gnsf);
 227
 228   /* Prepare reference frame */
 229   if (bPBC) {
 230       gpbc = gmx_rmpbc_init(&top->idef,ePBC,top->atoms.nr,box);
 231       gmx_rmpbc(gpbc,top->atoms.nr,box,x);
 232   }
 233
 234   natoms=read_first_x(oenv,&status,fnTRX,&t,&x,box);
 235   if (natoms != top->atoms.nr) {
 236       fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,top->atoms.nr);
 237   }
 238
 239   do {
 240       if (bPBC) {
 241           gmx_rmpbc(gpbc,top->atoms.nr,box,x);
 242       }
 243       /* allocate memory for pr */
 244       if (pr == NULL) {
 245           /* in case its first frame to read */
 246           snew(pr,1);
 247       }
 248       /*  realy calc p(r) */
 249       prframecurrent = calc_radial_distribution_histogram(gsans,x,box,index,isize,binwidth,bMC,bNORM,mcover,seed);
 250       /* copy prframecurrent -> pr and summ up pr->gr[i] */
 251       /* allocate and/or resize memory for pr->gr[i] and pr->r[i] */
 252       if (pr->gr == NULL) {
 253           /* check if we use pr->gr first time */
 254           snew(pr->gr,prframecurrent->grn);
 255           snew(pr->r,prframecurrent->grn);
 256       } else {
 257           /* resize pr->gr and pr->r if needed to preven overruns */
 258           if(prframecurrent->grn > pr->grn) {
 259               srenew(pr->gr,prframecurrent->grn);
 260               srenew(pr->r,prframecurrent->grn);
 261           }
 262       }
 263       pr->grn = prframecurrent->grn;
 264       pr->binwidth = prframecurrent->binwidth;
 265       /* summ up gr and fill r */
 266       for(i=0;i<prframecurrent->grn;i++) {
 267           pr->gr[i] += prframecurrent->gr[i];
 268           pr->r[i] = prframecurrent->r[i];
 269       }
 270       /* normalize histo */
 271       normalize_probability(prframecurrent->grn,prframecurrent->gr);
 272       /* convert p(r) to sq */
 273       sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent,start_q,end_q,q_step);
 274       /* print frame data if needed */
 275       if(opt2fn_null("-prframe",NFILE,fnm)) {
 276           snew(hdr,25);
 277           snew(suffix,GMX_PATH_MAX);
 278           /* prepare header */
 279           sprintf(hdr,"g(r), t = %f",t);
 280           /* prepare output filename */
 281           fnmdup = dup_tfn(NFILE,fnm);
 282           sprintf(suffix,"-t%.2f",t);
 283           add_suffix_to_output_names(fnmdup,NFILE,suffix);
 284           fp = xvgropen(opt2fn_null("-prframe",NFILE,fnmdup),hdr,"Distance (nm)","Probability",oenv);
 285           for(i=0;i<prframecurrent->grn;i++) {
 286               fprintf(fp,"%10.6f%10.6f\n",prframecurrent->r[i],prframecurrent->gr[i]);
 287           }
 288           done_filenms(NFILE,fnmdup);
 289           fclose(fp);
 290           sfree(hdr);
 291           sfree(suffix);
 292           sfree(fnmdup);
 293       }
 294       if(opt2fn_null("-sqframe",NFILE,fnm)) {
 295           snew(hdr,25);
 296           snew(suffix,GMX_PATH_MAX);
 297           /* prepare header */
 298           sprintf(hdr,"I(q), t = %f",t);
 299           /* prepare output filename */
 300           fnmdup = dup_tfn(NFILE,fnm);
 301           sprintf(suffix,"-t%.2f",t);
 302           add_suffix_to_output_names(fnmdup,NFILE,suffix);
 303           fp = xvgropen(opt2fn_null("-sqframe",NFILE,fnmdup),hdr,"q (nm^-1)","s(q)/s(0)",oenv);
 304           for(i=0;i<sqframecurrent->qn;i++) {
 305               fprintf(fp,"%10.6f%10.6f\n",sqframecurrent->q[i],sqframecurrent->s[i]);
 306           }
 307           done_filenms(NFILE,fnmdup);
 308           fclose(fp);
 309           sfree(hdr);
 310           sfree(suffix);
 311           sfree(fnmdup);
 312       }
 313       /* free pr structure */
 314       sfree(prframecurrent->gr);
 315       sfree(prframecurrent->r);
 316       sfree(prframecurrent);
 317       /* free sq structure */
 318       sfree(sqframecurrent->q);
 319       sfree(sqframecurrent->s);
 320       sfree(sqframecurrent);
 321   } while (read_next_x(oenv,status,&t,natoms,x,box));
 322   close_trj(status);
 323
 324   /* normalize histo */
 325   normalize_probability(pr->grn,pr->gr);
 326   sq = convert_histogram_to_intensity_curve(pr,start_q,end_q,q_step);
 327   /* prepare pr.xvg */
 328   fp = xvgropen(opt2fn_null("-pr",NFILE,fnm),"G(r)","Distance (nm)","Probability",oenv);
 329   for(i=0;i<pr->grn;i++)
 330       fprintf(fp,"%10.6f%10.6f\n",pr->r[i],pr->gr[i]);
 331   xvgrclose(fp);
 332
 333   /* prepare sq.xvg */
 334   fp = xvgropen(opt2fn_null("-sq",NFILE,fnm),"I(q)","q (nm^-1)","s(q)/s(0)",oenv);
 335   for(i=0;i<sq->qn;i++) {
 336       fprintf(fp,"%10.6f%10.6f\n",sq->q[i],sq->s[i]);
 337   }
 338   xvgrclose(fp);
 339   /*
 340    * Clean up memory
 341    */
 342   sfree(pr->gr);
 343   sfree(pr->r);
 344   sfree(pr);
 345   sfree(sq->q);
 346   sfree(sq->s);
 347   sfree(sq);
 348
 349   please_cite(stdout,"Garmay2012");
 350   thanx(stderr);
 351
 352   return 0;
 353 }