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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32 static char *SRCID_disco_c = "$Id$";
  33 #include "macros.h"
  34 #include "statutil.h"
  35 #include "pdbio.h"
  36 #include "smalloc.h"
  37 #include "random.h"
  38 #include "vec.h"
  39 #include "princ.h"
  40 #include "confio.h"
  41 #include "rdgroup.h"
  42 #include "filenm.h"
  43 #include "do_fit.h"
  44 #include "tpxio.h"
  45 #include "copyrite.h"
  46 #include "disco.h"
  47 #include "xvgr.h"
  48 #include "main.h"
  49 #include "network.h"
  50
  51 void rand_box(bool bUserBox,
  52               matrix box,rvec boxsize,int nres,bool bCubic,int *seed)
  53 {
  54   int  m;
  55   real fac;
  56
  57   clear_mat(box);
  58
  59   if (bUserBox) {
  60     for(m=0; (m<DIM); m++)
  61       box[m][m] = boxsize[m];
  62   }
  63   else {
  64     /* Generate a random box with size between 5*nres and 10*nres in nm */
  65     fac = 0.5*nres; /* Ca-Ca distance is 0.35 nm */
  66     box[XX][XX] = fac*(1+rando(seed));
  67     if (bCubic)
  68       box[YY][YY] = box[ZZ][ZZ] = box[XX][XX];
  69     else {
  70       box[YY][YY] = fac*(1+rando(seed));
  71       box[ZZ][ZZ] = fac*(1+rando(seed));
  72     }
  73     for(m=0; (m<DIM); m++)
  74       boxsize[m] = box[m][m];
  75   }
  76 }
  77
  78 void rand_coord(rvec x,int *seed,rvec box)
  79 {
  80   int m;
  81
  82   for(m=0; (m<DIM); m++)
  83     x[m] = box[m]*rando(seed);
  84 }
  85
  86 void rand_coords(int natom,rvec x[],rvec xref[],real weight[],
  87                  bool bCenter,rvec xcenter[],rvec box,int *seed)
  88 {
  89   int i;
  90
  91   for(i=0; (i<natom); i++) {
  92     if (weight[i] == 0)
  93       copy_rvec(xref[i],x[i]);
  94     else {
  95       rand_coord(x[i],seed,box);
  96       if (bCenter)
  97         rvec_inc(x[i],xcenter[i]);
  98     }
  99   }
 100 }
 101
 102 void pr_conv_stat(FILE *fp,int ntry,int nconv,double tnit)
 103 {
 104   fprintf(fp,"\n-------------------------\n");
 105   fprintf(fp," Convergence statistics:\n");
 106   fprintf(fp," # tries:     %d\n",ntry);
 107   fprintf(fp," # converged: %d\n",nconv);
 108   fprintf(fp," # nit/ntry:  %d\n",(int)(tnit/ntry));
 109   if (nconv > 0)
 110     fprintf(fp," # nit/nconv: %d\n",(int)(tnit/nconv));
 111   fprintf(fp,"-------------------------\n");
 112 }
 113
 114 static void do_disco(FILE *log,char *outfn,char *keepfn,t_correct *c,
 115                      bool bVerbose,t_atoms *atoms,
 116                      rvec xref[],rvec xcenter[],
 117                      int nstruct,int *seed,
 118                      bool bFit,int nfit,atom_id fit_ind[],
 119                      bool bPrintViol,char *violfn,rvec boxsize)
 120 {
 121   FILE    *fp,*gp;
 122   int     *nconvdist;
 123   int     i,k,kk,nconv,ntry,status,kstatus,natom,nres,nit,nvtest,nviol;
 124   double  tnit;
 125   rvec    *x,xcm;
 126   matrix  box,wrbox;
 127   atom_id *wr_ind;
 128   real    *w_rls;
 129   bool    bConverged;
 130
 131   natom  = atoms->nr;
 132   nres   = atoms->nres;
 133
 134   /* Initiate some local arrays */
 135   init_corr2(c,natom);
 136
 137   clear_mat(wrbox);
 138   wrbox[XX][XX] = wrbox[YY][YY] = wrbox[ZZ][ZZ] = nres;
 139   status = open_trx(outfn,"w");
 140   if (keepfn)
 141     kstatus = open_trx(keepfn,"w");
 142   else
 143     kstatus = -1;
 144   snew(x,natom);
 145   snew(wr_ind,natom);
 146   for(k=0; (k<natom); k++)
 147     wr_ind[k]=k;
 148   snew(w_rls,natom);
 149   for(k=0; (k<nfit); k++)
 150     w_rls[fit_ind[k]] = 1;
 151
 152   snew(nconvdist,c->maxnit+1);
 153   /* Now loop over structures */
 154   tnit = 0;
 155   for(k=nconv=ntry=0; (k<nstruct); ntry++) {
 156     if (bVerbose)
 157       fprintf(stderr,"\rTry: %d, Success: %d",ntry,nconv);
 158
 159     /* Generate random box*/
 160     rand_box(c->bBox,box,boxsize,nres,c->bCubic,seed);
 161
 162     /* Generate random coords */
 163     rand_coords(natom,x,xref,c->weight,c->bCenter,xcenter,boxsize,seed);
 164
 165     /* Now correct the random coords */
 166     nviol = shake_coords(log,bVerbose,k,natom,xref,x,seed,box,c,&nit);
 167     bConverged = (nviol == 0);
 168     tnit += nit;
 169
 170     if (bConverged)
 171       nconvdist[nit]++;
 172
 173     nvtest = quick_check(bVerbose ? log : NULL,natom,x,box,c);
 174     fprintf(stderr,"Double checking: %d violations\n",nvtest);
 175
 176     if (bConverged || keepfn) {
 177       center_in_box(natom,x,wrbox,x);
 178       if (bFit)
 179         do_fit(natom,w_rls,xref,x);
 180       write_trx(bConverged ? status : kstatus,
 181                 natom,wr_ind,atoms,k,(real) k,wrbox,x,NULL);
 182
 183       if (bConverged)
 184         nconv++;
 185
 186       k++;
 187     }
 188     if (bPrintViol) {
 189       /* Print structure coloured by the violations */
 190       if (!atoms->pdbinfo)
 191         snew(atoms->pdbinfo,natom);
 192       for(kk=0; (kk<natom); kk++)
 193         atoms->pdbinfo[kk].bfac = (real) c->bViol[kk];
 194       gp=ffopen(violfn,"w");
 195       write_pdbfile(gp,"Structure coloured by violation",atoms,x,box,0,-1);
 196       ffclose(gp);
 197     }
 198   }
 199   close_trx(status);
 200   if (keepfn)
 201     close_trx(kstatus);
 202   gp = xvgropen("conv_stat.xvg","Iterations per converged structure",
 203                 "nit","N");
 204   for(i=0; (i<c->maxnit); i++)
 205     fprintf(gp,"%10d  %10d\n",i,nconvdist[i]);
 206   ffclose(gp);
 207   sfree(x);
 208   sfree(w_rls);
 209   sfree(wr_ind);
 210   sfree(nconvdist);
 211
 212   pr_conv_stat(log,ntry,nconv,tnit);
 213   pr_conv_stat(stderr,ntry,nconv,tnit);
 214 }
 215
 216 int main(int argc,char *argv[])
 217 {
 218   static char *desc[] = {
 219     "disco reads a topology (tpr) file and runs distance geometry",
 220     "calculations based on the distances defined in the",
 221     "distance-restraints section of the topology. An appropriate tpr",
 222     "file may be generated by the cdist program.[PAR]",
 223     "The algorithm is the CONCOORD algorithm of De Groot et al.,",
 224     "which in turn is derived from the SHAKE alogrithm.[PAR]",
 225     "A parallel version of disco is under development whihc uses a",
 226     "master-slave approach. Slaves work asynchronously, and it is no",
 227     "problem when nodes are not equally fast, or when a node dies,",
 228     "unless it is the master node."
 229   };
 230   FILE        *fp,*dp;
 231   char        title[256];
 232   bool        bCenter,bBox;
 233   t_commrec   *cr;
 234   t_atoms     atoms,newatoms;
 235   t_correct   *corr;
 236   rvec        xcm,*xref,*xcenter=NULL;
 237   matrix      box;
 238   real        t,lambda,tot_weight;
 239   int         i,nfit,step,natom;
 240   atom_id     *fit_ind;
 241   char        *grpname;
 242
 243   static int  nstruct=10,maxnit=1000,seed=1997,nbcheck=1;
 244   static int  nstprint=1,nstranlist=0,ngrow=0;
 245   static bool bVerbose=TRUE,bCubic=FALSE,bWeight=FALSE,bLower=FALSE;
 246   static real lowdev=0.05;
 247   static bool bExplicit=FALSE,bChiral=TRUE,bFit=FALSE,bDump=FALSE,bPep=TRUE;
 248   static bool bRanlistFirst=TRUE;
 249   static rvec boxsize={ 2, 2, 2 };
 250   t_pargs pa[] = {
 251     { "-nf",    FALSE, etINT,     {&nstruct},
 252       "Number of structures to generate" },
 253     { "-nit",   FALSE, etINT,     {&maxnit},
 254       "Max number of iterations for a structure to converge" },
 255     { "-v",     FALSE, etBOOL,    {&bVerbose},
 256       "Be verbosive" },
 257     { "-chiral",   FALSE, etBOOL, {&bChiral},
 258       "Check chirality during disco-ing" },
 259     { "-pep",   FALSE,  etBOOL,   {&bPep},
 260       "Flip all cis-peptide bonds automatically to trans" },
 261     { "-lower", FALSE,  etBOOL,   {&bLower},
 262       "Use lower bounds only for nonbondeds." },
 263     { "-weighted", FALSE, etBOOL, {&bWeight},
 264       "Use weighted disco. The STX file must be a pdb file in this case and weights are read from the occupancy field" },
 265     { "-dump",     FALSE, etBOOL, {&bDump},
 266       "Dump the trajectory of the shaking to testX.xtc file where X is the structure number." },
 267     { "-cubic",    FALSE, etBOOL, {&bCubic},
 268       "Generate coordinates in a cubic box, rather than rectangular" },
 269     { "-explicit", FALSE, etBOOL, {&bExplicit},
 270       "Use explicit updating of positions if the sum of deviations is smaller than lowdev" },
 271     { "-fit",      FALSE, etBOOL, {&bFit},
 272       "Fit output structures to reference structure in tpx file" },
 273     { "-nbcheck",  FALSE, etINT,  {&nbcheck},
 274       "Check non-bonded interactions every N steps" },
 275     { "-nstprint", FALSE, etINT,  {&nstprint},
 276       "Print number of violations every N steps" },
 277     { "-ranlist",  FALSE, etINT,  {&nstranlist},
 278       "Update list order to avoid bias every n steps" },
 279     { "-ranlistfirst",  FALSE, etBOOL,  {&bRanlistFirst},
 280       "Randomize list once before shaking" },
 281     { "-lowdev",   FALSE, etREAL, {&lowdev},
 282       "Low deviation [Sum of distance deviation per atom in nm] beyond which nonbondeds are done every step" },
 283     { "-seed",     FALSE, etINT,  {&seed},
 284       "Seed for the random number generator" },
 285     { "-box",      FALSE, etRVEC, {boxsize},
 286       "Boxsize (nm) for generating random coordinates" },
 287     { "-grow",     FALSE, etINT,  {&ngrow},
 288       "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" }
 289
 290   };
 291 #define NPA asize(pa)
 292
 293   t_filenm fnm[] = {
 294     { efLOG, "-g",     "disco", ffWRITE },
 295     { efSTX, "-f",      NULL,   ffREAD  },
 296     { efDAT, "-d",     "cdist", ffREAD },
 297     { efDAT, "-do",    "distout",  ffOPTWR },
 298     { efSTO, "-c",      NULL,   ffREAD  },
 299     { efSTO, "-center", NULL,   ffOPTRD },
 300     { efNDX, "-n",      NULL,   ffOPTRD },
 301     { efTRX, "-o",   "structs", ffWRITE },
 302     { efTRX, "-keep", "unconverged", ffOPTWR },
 303     { efPDB, "-viol",  "vvv",   ffOPTWR }
 304   };
 305 #define NFILE asize(fnm)
 306
 307   cr       = init_par(&argc,&argv);
 308   bVerbose = bVerbose && MASTER(cr);
 309
 310   if (MASTER(cr)) {
 311     CopyRight(stderr,argv[0]);
 312
 313     parse_common_args(&argc,argv,PCA_BE_NICE,
 314                       NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 315
 316     /* Copy arguments to correct structure */
 317     bCenter = opt2bSet("-center",NFILE,fnm);
 318     bBox    = opt2parg_bSet("-box",asize(pa),pa),
 319     corr    = init_corr(maxnit,nstprint,nbcheck,nstranlist,ngrow,bExplicit,
 320                         bChiral,bPep,bDump,lowdev,bLower,bRanlistFirst,
 321                         bCubic,bBox,bCenter);
 322   }
 323   /* Open the log file */
 324   open_log(ftp2fn(efLOG,NFILE,fnm),cr);
 325
 326   if (MASTER(cr)) {
 327     CopyRight(stdlog,argv[0]);
 328     please_cite(stdlog,"Ryckaert77a");
 329     please_cite(stdlog,"DeGroot97a");
 330   }
 331
 332   if (MASTER(cr)) {
 333     /* Get number of atoms etc. */
 334     get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natom);
 335
 336     init_t_atoms(&atoms,natom,bWeight);
 337     snew(xref,natom);
 338     read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,xref,NULL,box);
 339
 340     /* Translate reference and xcenter coords to C.O.M. */
 341     sub_xcm(xref,natom,NULL,NULL,xcm,FALSE);
 342
 343     snew(corr->weight,natom);
 344     tot_weight = 0;
 345     for(i=0; (i<natom); i++) {
 346       corr->weight[i] = bWeight ? atoms.pdbinfo[i].occup : 1;
 347       tot_weight += corr->weight[i];
 348     }
 349
 350     fprintf(stderr,"Reading distances from %s\n",opt2fn("-d",NFILE,fnm));
 351     read_dist(stdlog,opt2fn("-d",NFILE,fnm),natom,corr);
 352
 353     /* Dump a distance file if necessary */
 354     if (opt2bSet("-do",NFILE,fnm)) {
 355       dp = fopen(opt2fn("-do",NFILE,fnm),"w");
 356       pr_distances(dp,corr);
 357       fclose(dp);
 358     }
 359
 360     /* Check distances */
 361     check_dist(stdlog,corr);
 362
 363     /* Read index if necessary */
 364     if (bFit) {
 365       fprintf(stderr,"Select group for fitting output structures:\n");
 366       get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,
 367                 &nfit,&fit_ind,&grpname);
 368     }
 369     else {
 370       nfit = 0;
 371       fit_ind = NULL;
 372     }
 373
 374     /* Read centers for generating coordinates (optional) */
 375     if (bCenter) {
 376       snew(xcenter,natom);
 377       init_t_atoms(&newatoms,natom,TRUE);
 378       read_stx_conf(opt2fn("-center",NFILE,fnm),title,
 379                     &newatoms,xcenter,NULL,box);
 380       free_t_atoms(&newatoms);
 381
 382       for(i=0; (i<natom); i++)
 383         rvec_dec(xcenter[i],xcm);
 384     }
 385
 386     /*
 387      * define improper dihedrals that are not automatically correct
 388      * when all distances are correct
 389      */
 390     define_impropers(stdlog,&atoms,corr);
 391
 392     /* define peptide-bonds, so we can correct cis to trans
 393      * Adam Kirrander 990121
 394      */
 395     define_peptide_bonds(stdlog,&atoms,corr);
 396
 397     /* Print parameters */
 398     pr_corr(stdlog,corr);
 399   }
 400
 401   /* Now do my thing */
 402 #ifdef USE_MPI
 403   if (PAR(cr)) {
 404     if (MASTER(cr))
 405       disco_master(cr,stdlog,opt2fn("-o",NFILE,fnm),
 406                    opt2bSet("-keep",NFILE,fnm) ? opt2fn("-keep",NFILE,fnm) : NULL,
 407                    corr,bVerbose,&atoms,
 408                    xref,xcenter,nstruct,&seed,bFit,nfit,fit_ind,
 409                    opt2bSet("-viol",NFILE,fnm),opt2fn("-viol",NFILE,fnm),boxsize);
 410     else
 411       disco_slave(cr,stdlog);
 412   }
 413   else {
 414 #endif
 415     do_disco(stdlog,opt2fn("-o",NFILE,fnm),
 416              opt2bSet("-keep",NFILE,fnm) ? opt2fn("-keep",NFILE,fnm) : NULL,
 417              corr,bVerbose,&atoms,
 418              xref,xcenter,nstruct,&seed,bFit,nfit,fit_ind,
 419              opt2bSet("-viol",NFILE,fnm),opt2fn("-viol",NFILE,fnm),
 420              boxsize);
 421 #ifdef USE_MPI
 422   }
 423 #endif
 424   ffclose(stdlog);
 425 #ifdef USE_MPI
 426   if (PAR(cr))
 427     MPI_Finalize();
 428 #endif
 429   if (MASTER(cr))
 430     thanx(stderr);
 431
 432   return 0;
 433 }