src/tools/disco.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_disco_c = "$Id$";
  30
  31 #include "macros.h"
  32 #include "statutil.h"
  33 #include "pdbio.h"
  34 #include "smalloc.h"
  35 #include "random.h"
  36 #include "vec.h"
  37 #include "princ.h"
  38 #include "confio.h"
  39 #include "rdgroup.h"
  40 #include "filenm.h"
  41 #include "do_fit.h"
  42 #include "tpxio.h"
  43 #include "copyrite.h"
  44 #include "disco.h"
  45 #include "xvgr.h"
  46
  47 void rand_box(bool bUserBox,
  48               matrix box,rvec boxsize,int nres,bool bCubic,int *seed)
  49 {
  50   int  m;
  51   real fac;
  52
  53   clear_mat(box);
  54
  55   if (bUserBox) {
  56     for(m=0; (m<DIM); m++)
  57       box[m][m] = boxsize[m];
  58   }
  59   else {
  60     /* Generate a random box with size between 5*nres and 10*nres in nm */
  61     fac = 0.5*nres; /* Ca-Ca distance is 0.35 nm */
  62     box[XX][XX] = fac*(1+rando(seed));
  63     if (bCubic)
  64       box[YY][YY] = box[ZZ][ZZ] = box[XX][XX];
  65     else {
  66       box[YY][YY] = fac*(1+rando(seed));
  67       box[ZZ][ZZ] = fac*(1+rando(seed));
  68     }
  69     for(m=0; (m<DIM); m++)
  70       boxsize[m] = box[m][m];
  71   }
  72 }
  73
  74 void rand_coord(rvec x,int *seed,rvec box)
  75 {
  76   int m;
  77
  78   for(m=0; (m<DIM); m++)
  79     x[m] = box[m]*rando(seed);
  80 }
  81
  82 void rand_coords(int natom,rvec x[],rvec xref[],real weight[],
  83                  bool bCenter,rvec xcenter[],rvec box,int *seed)
  84 {
  85   int i;
  86
  87   for(i=0; (i<natom); i++) {
  88     if (weight[i] == 0)
  89       copy_rvec(xref[i],x[i]);
  90     else {
  91       rand_coord(x[i],seed,box);
  92       if (bCenter)
  93         rvec_inc(x[i],xcenter[i]);
  94     }
  95   }
  96 }
  97
  98 static void pr_conv_stat(FILE *fp,int ntry,int nconv,double tnit)
  99 {
 100   fprintf(fp,"\n-------------------------\n");
 101   fprintf(fp," Convergence statistics:\n");
 102   fprintf(fp," # tries:     %d\n",ntry);
 103   fprintf(fp," # converged: %d\n",nconv);
 104   fprintf(fp," # nit/ntry:  %d\n",(int)(tnit/ntry));
 105   if (nconv > 0)
 106     fprintf(fp," # nit/nconv: %d\n",(int)(tnit/nconv));
 107   fprintf(fp,"-------------------------\n");
 108 }
 109
 110
 111 void do_disco(FILE *log,char *outfn,char *keepfn,t_correct *c,
 112               bool bVerbose,t_atoms *atoms,
 113               rvec xref[],bool bCenter,rvec xcenter[],real weight[],
 114               int nstruct,bool bCubic,int *seed,
 115               bool bFit,int nfit,atom_id fit_ind[],
 116               bool bPrintViol,char *violfn,
 117               bool bBox,rvec boxsize)
 118 {
 119   FILE    *fp,*gp;
 120   int     *nconvdist;
 121   int     i,k,kk,nconv,ntry,status,kstatus,natom,nres,nit,nvtest;
 122   double  tnit;
 123   rvec    *x,xcm;
 124   matrix  box,wrbox;
 125   atom_id *wr_ind;
 126   real    *w_rls;
 127   bool    bConverged;
 128
 129   natom  = atoms->nr;
 130   nres   = atoms->nres;
 131   wrbox[XX][XX] = wrbox[YY][YY] = wrbox[ZZ][ZZ] = nres;
 132   status = open_trx(outfn,"w");
 133   if (keepfn)
 134     kstatus = open_trx(keepfn,"w");
 135   else
 136     kstatus = -1;
 137   snew(x,natom);
 138   snew(wr_ind,natom);
 139   for(k=0; (k<natom); k++)
 140     wr_ind[k]=k;
 141   snew(w_rls,natom);
 142   for(k=0; (k<nfit); k++)
 143     w_rls[fit_ind[k]] = 1;
 144
 145   snew(nconvdist,c->maxnit+1);
 146   /* Now loop over structures */
 147   tnit = 0;
 148   for(k=nconv=ntry=0; (k<nstruct); ntry++) {
 149     if (bVerbose)
 150       fprintf(stderr,"\rTry: %d, Success: %d",ntry,nconv);
 151
 152     /* Generate random box*/
 153     rand_box(bBox,box,boxsize,nres,bCubic,seed);
 154
 155     /* Generate random coords */
 156     rand_coords(natom,x,xref,weight,bCenter,xcenter,boxsize,seed);
 157
 158     /* Now correct the random coords */
 159     bConverged = shake_coords(log,bVerbose,k,natom,xref,x,seed,box,c,&nit);
 160     tnit += nit;
 161
 162     if (bConverged)
 163       nconvdist[nit]++;
 164
 165     nvtest = quick_check(bVerbose ? log : NULL,natom,x,box,c);
 166     fprintf(stderr,"Double checking: %d violations\n",nvtest);
 167
 168     if (bConverged || keepfn) {
 169       center_in_box(natom,x,wrbox,x);
 170       if (bFit)
 171         do_fit(natom,w_rls,xref,x);
 172       write_trx(bConverged ? status : kstatus,
 173                 natom,wr_ind,atoms,k,(real) k,wrbox,x,NULL);
 174
 175       if (bConverged)
 176         nconv++;
 177
 178       k++;
 179     }
 180     if (bPrintViol) {
 181       /* Print structure coloured by the violations */
 182       if (!atoms->pdbinfo)
 183         snew(atoms->pdbinfo,natom);
 184       for(kk=0; (kk<natom); kk++)
 185         atoms->pdbinfo[kk].bfac = (real) c->bViol[kk];
 186       gp=ffopen(violfn,"w");
 187       write_pdbfile(gp,"Structure coloured by violation",
 188                     atoms,x,box,'A',TRUE);
 189       ffclose(gp);
 190     }
 191   }
 192   close_trx(status);
 193   if (keepfn)
 194     close_trx(kstatus);
 195   gp = xvgropen("conv_stat.xvg","Iterations per converged structure",
 196                 "nit","N");
 197   for(i=0; (i<c->maxnit); i++)
 198     fprintf(gp,"%10d  %10d\n",i,nconvdist[i]);
 199   ffclose(gp);
 200   sfree(x);
 201   sfree(w_rls);
 202   sfree(wr_ind);
 203   sfree(nconvdist);
 204
 205   pr_conv_stat(log,ntry,nconv,tnit);
 206   pr_conv_stat(stderr,ntry,nconv,tnit);
 207 }
 208
 209 int main(int argc,char *argv[])
 210 {
 211   static char *desc[] = {
 212     "disco reads a topology (tpr) file and runs distance geometry",
 213     "calculations based on the distances defined in the",
 214     "distance-restraints section of the topology. An appropriate tpr",
 215     "file may be generated by the cdist program.[PAR]",
 216     "The algorithm is the CONCOORD algorithm of De Groot et al.,",
 217     "which in turn is derived from the SHAKE alogrithm"
 218   };
 219   FILE        *fp,*dp;
 220   char        title[256];
 221   bool        bCenter;
 222   t_atoms     atoms,newatoms;
 223   t_correct   *corr;
 224   rvec        xcm,*xref,*xcenter=NULL;
 225   matrix      box;
 226   real        *weight,t,lambda,tot_weight;
 227   int         i,nfit,step,natom;
 228   atom_id     *fit_ind;
 229   char        *grpname;
 230
 231   static int  nstruct=10,maxnit=1000,seed=1997,nbcheck=1;
 232   static int  nstprint=1,nstranlist=1,ngrow=0;
 233   static bool bVerbose=TRUE,bCubic=FALSE,bWeight=FALSE,bLower=FALSE;
 234   static real lowdev=0.05,cutoff=0;
 235   static bool bExplicit=FALSE,bChiral=TRUE,bFit=FALSE,bDump=FALSE,bPep=TRUE;
 236   static rvec boxsize={ 2, 2, 2 };
 237   t_pargs pa[] = {
 238     { "-nf",    FALSE, etINT,     &nstruct,
 239       "Number of structures to generate" },
 240     { "-nit",   FALSE, etINT,     &maxnit,
 241       "Max number of iterations for a structure to converge" },
 242     { "-v",     FALSE, etBOOL,    &bVerbose,
 243       "Be verbosive" },
 244     { "-chiral",   FALSE, etBOOL, &bChiral,
 245       "Check chirality during disco-ing" },
 246     { "-pep",   FALSE,  etBOOL,   &bPep,
 247       "Flip all cis-peptide bonds automatically to trans" },
 248     { "-lower", FALSE,  etBOOL,   &bLower,
 249       "Use lower bounds only for nonbondeds." },
 250     { "-weighted", FALSE, etBOOL, &bWeight,
 251       "Use weighted disco. The STX file must be a pdb file in this case and weights are read from the occupancy field" },
 252     { "-cutoff",   FALSE, etREAL, &cutoff,
 253       "Cut-off for taking pairs into account when measuring distance" },
 254     { "-dump",     FALSE, etBOOL, &bDump,
 255       "Dump the trajectory of the shaking to testX.xtc file where X is the structure number." },
 256     { "-cubic",    FALSE, etBOOL, &bCubic,
 257       "Generate coordinates in a cubic box, rather than rectangular" },
 258     { "-explicit", FALSE, etBOOL, &bExplicit,
 259       "Use explicit updating of positions if the sum of deviations is smaller than lowdev" },
 260     { "-fit",      FALSE, etBOOL, &bFit,
 261       "Fit output structures to reference structure in tpx file" },
 262     { "-nbcheck",  FALSE, etINT,  &nbcheck,
 263       "Check non-bonded interactions every N steps" },
 264     { "-nstprint", FALSE, etINT,  &nstprint,
 265       "Print number of violations every N steps" },
 266     { "-ranlist",  FALSE, etINT,  &nstranlist,
 267       "Update list order to avoid bias every n steps" },
 268     { "-lowdev",   FALSE, etREAL, &lowdev,
 269       "Low deviation [Sum of distance deviation per atom in nm] beyond which nonbondeds are done every step" },
 270     { "-seed",     FALSE, etINT,  &seed,
 271       "Seed for the random number generator" },
 272     { "-box",      FALSE, etRVEC, boxsize,
 273       "Boxsize (nm) for generating random coordinates" },
 274     { "-grow",     FALSE, etINT,  &ngrow,
 275       "Number of steps after which Van der Waals lower bounds grow from 0 to the real lower bounds. If this is 0 (default), the Van der Waals lower bounds are in effect from the beginning" }
 276
 277   };
 278 #define NPA asize(pa)
 279
 280   t_filenm fnm[] = {
 281     { efLOG, "-g",     "disco", ffWRITE },
 282     { efSTX, "-f",      NULL,   ffREAD  },
 283     { efDAT, "-d",     "cdist", ffREAD },
 284     { efDAT, "-do",    "distout",  ffOPTWR },
 285     { efSTO, "-c",      NULL,   ffREAD  },
 286     { efSTO, "-center", NULL,   ffOPTRD },
 287     { efNDX, "-n",      NULL,   ffOPTRD },
 288     { efTRX, "-o",   "structs", ffWRITE },
 289     { efTRX, "-keep", "unconverged", ffOPTWR },
 290     { efPDB, "-viol",  "vvv",   ffOPTWR }
 291   };
 292 #define NFILE asize(fnm)
 293
 294   CopyRight(stdout,argv[0]);
 295
 296   parse_common_args(&argc,argv,0,TRUE,
 297                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 298   /* Copy arguments to correct structure */
 299   corr = init_corr(maxnit,nstprint,nbcheck,nstranlist,ngrow,bExplicit,
 300                    bChiral,bPep,bDump,lowdev,bLower);
 301
 302   /* Open the log file */
 303   fp = ftp2FILE(efLOG,NFILE,fnm,"w");
 304
 305   please_cite(fp,"Ryckaert77a");
 306   please_cite(fp,"DeGroot97a");
 307
 308   init_lookup_table(fp);
 309
 310   /* Get number of atoms etc. */
 311   get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natom);
 312
 313   init_t_atoms(&atoms,natom,bWeight);
 314   snew(xref,natom);
 315   read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms, xref,NULL,box);
 316
 317   snew(weight,natom);
 318   tot_weight = 0;
 319   for(i=0; (i<natom); i++) {
 320     weight[i] = bWeight ? atoms.pdbinfo[i].occup : 1;
 321     tot_weight += weight[i];
 322   }
 323
 324   fprintf(stderr,"Reading distances from %s\n",opt2fn("-d",NFILE,fnm));
 325   read_dist(fp,opt2fn("-d",NFILE,fnm),natom,corr,weight);
 326
 327   /* Dump a distance file if necessary */
 328   if (opt2bSet("-do",NFILE,fnm)) {
 329     dp = fopen(opt2fn("-do",NFILE,fnm),"w");
 330     pr_distances(dp,corr);
 331     fclose(dp);
 332   }
 333
 334   /* Check distances */
 335   check_dist(fp,corr);
 336
 337   /* Make tags */
 338   make_tags(corr,natom);
 339
 340   /* Translate reference and xcenter coords to C.O.M. */
 341   sub_xcm(xref,natom,NULL,NULL,xcm,FALSE);
 342
 343   /* Read index if necessary */
 344   if (bFit) {
 345     fprintf(stderr,"Select group for fitting output structures:\n");
 346     get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 347               &nfit,&fit_ind,&grpname);
 348   }
 349   else {
 350     nfit = 0;
 351     fit_ind = NULL;
 352   }
 353
 354   /* Read centers for generating coordinates (optional) */
 355   bCenter = opt2bSet("-center",NFILE,fnm);
 356   if (bCenter) {
 357     snew(xcenter,natom);
 358     init_t_atoms(&newatoms,natom,TRUE);
 359     read_stx_conf(opt2fn("-center",NFILE,fnm),title,
 360                   &newatoms,xcenter,NULL,box);
 361     free_t_atoms(&newatoms);
 362
 363     for(i=0; (i<natom); i++)
 364       rvec_dec(xcenter[i],xcm);
 365   }
 366
 367   /*
 368    * define improper dihedrals that are not automatically correct
 369    * when all distances are correct
 370    */
 371   define_impropers(fp,&atoms,corr);
 372
 373   /* define peptide-bonds, so we can correct cis to trans
 374    * Adam Kirrander 990121
 375    */
 376   define_peptide_bonds(fp,&atoms,corr);
 377
 378   /* Print parameters */
 379   pr_corr(fp,corr);
 380
 381   /* Now do my thing */
 382   do_disco(fp,opt2fn("-o",NFILE,fnm),
 383            opt2bSet("-keep",NFILE,fnm) ? opt2fn("-keep",NFILE,fnm) : NULL,
 384            corr,bVerbose,&atoms,
 385            xref,bCenter,xcenter,weight,nstruct,bCubic,&seed,bFit,nfit,fit_ind,
 386            opt2bSet("-viol",NFILE,fnm),opt2fn("-viol",NFILE,fnm),
 387            opt2parg_bSet("-box",asize(pa),pa),boxsize);
 388   ffclose(fp);
 389
 390   thanx(stdout);
 391
 392   return 0;
 393 }