src/tools/discopar.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
  33 #include "typedefs.h"
  34 #include "network.h"
  35 #include "smalloc.h"
  36 #include "vec.h"
  37 #include "statutil.h"
  38 #include "do_fit.h"
  39 #include "random.h"
  40 #include "futil.h"
  41 #include "xvgr.h"
  42 #include "pdbio.h"
  43 #include "disco.h"
  44
  45 static t_correct *recv_init(FILE *fp,
  46                             t_commrec *cr,int *seed,int *natom,int *nres,
  47                             rvec **xref,rvec **xcenter,bool *bKeep)
  48 {
  49   t_correct *c;
  50
  51   gmx_rxs(0,seed,sizeof(*seed));
  52   gmx_rxs(0,natom,sizeof(*natom));
  53   gmx_rxs(0,nres,sizeof(*nres));
  54   gmx_rxs(0,bKeep,sizeof(*bKeep));
  55
  56   snew(c,1);
  57 #define cget(nn) gmx_rxs(0,record(c->nn))
  58   cget(maxnit);
  59   cget(nbcheck);
  60   cget(nstprint);
  61   cget(nstranlist);
  62   cget(ngrow);
  63   cget(bExplicit);
  64   cget(bChiral);
  65   cget(bPep);
  66   cget(bLowerOnly);
  67   cget(bRanlistFirst);
  68   cget(bCubic);
  69   cget(bBox);
  70   cget(bCenter);
  71   cget(lodev);
  72   cget(maxdist);
  73   cget(ndist);
  74   cget(npep);
  75   cget(nimp);
  76 #undef cget
  77   snew(*xref,*natom);
  78   gmx_rxs(0,*xref,*natom*sizeof(**xref));
  79   if (c->bCenter) {
  80     snew(*xcenter,*natom);
  81     gmx_rxs(0,*xcenter,*natom*sizeof(**xcenter));
  82   }
  83
  84   /* Get the important input data */
  85   snew(c->d,c->ndist);
  86   gmx_rxs(0,array(c->d,c->ndist));
  87
  88   snew(c->pepbond,c->npep);
  89   gmx_rxs(0,array(c->pepbond,c->npep));
  90
  91   snew(c->imp,c->nimp+1);
  92   gmx_rxs(0,array(c->imp,c->nimp+1));
  93
  94   snew(c->weight,*natom);
  95   gmx_rxs(0,array(c->weight,*natom));
  96
  97   /* Other arrays can be deduced from these */
  98   fprintf(fp,"Succesfully got input data\n");
  99   fflush(fp);
 100
 101   return c;
 102 }
 103
 104 static void send_init(FILE *fp,
 105                       t_commrec *cr,t_correct *c,int *seed,int natom,int nres,
 106                       rvec *xref,rvec *xcenter,bool bKeep)
 107 {
 108   int nodeid;
 109
 110   for(nodeid=1; (nodeid < cr->nnodes); nodeid++) {
 111     gmx_txs(nodeid,seed,sizeof(*seed));
 112     /* Update random seed */
 113     (void) rando(seed);
 114     gmx_txs(nodeid,record(natom));
 115     gmx_txs(nodeid,record(nres));
 116     gmx_txs(nodeid,record(bKeep));
 117
 118 #define cput(nn) gmx_txs(nodeid,record(c->nn))
 119     cput(maxnit);
 120     cput(nbcheck);
 121     cput(nstprint);
 122     cput(nstranlist);
 123     cput(ngrow);
 124     cput(bExplicit);
 125     cput(bChiral);
 126     cput(bPep);
 127     cput(bLowerOnly);
 128     cput(bRanlistFirst);
 129     cput(bCubic);
 130     cput(bBox);
 131     cput(bCenter);
 132     cput(lodev);
 133     cput(maxdist);
 134     cput(ndist);
 135     cput(npep);
 136     cput(nimp);
 137 #undef cput
 138     gmx_txs(nodeid,array(xref,natom));
 139
 140     if (c->bCenter)
 141       gmx_txs(nodeid,array(xcenter,natom));
 142
 143     /* Send the important input data */
 144     gmx_txs(nodeid,array(c->d,c->ndist));
 145
 146     gmx_txs(nodeid,array(c->pepbond,c->npep));
 147
 148     gmx_txs(nodeid,array(c->imp,c->nimp+1));
 149
 150     gmx_txs(nodeid,array(c->weight,natom));
 151
 152     /* Other arrays can be deduced from these */
 153     fprintf(fp,"Succesfully sent input data to nodeid %d\n",nodeid);
 154     fflush(fp);
 155   }
 156 }
 157
 158 static bool send_coords(t_commrec *cr,int nviol,int nit,int k,
 159                         int natom,rvec x[],matrix box,bool bKeep)
 160 {
 161   bool bDone;
 162
 163   debug_gmx();
 164   gmx_tx(0,record(cr->nodeid));
 165   gmx_tx_wait(0);
 166
 167   gmx_rxs(0,record(bDone));
 168   if (!bDone) {
 169     debug_gmx();
 170     gmx_txs(0,record(nviol));
 171     gmx_txs(0,record(nit));
 172     gmx_txs(0,record(k));
 173     if (bKeep || (nviol == 0)) {
 174       gmx_txs(0,record(natom));
 175       gmx_txs(0,array(x,natom));
 176       gmx_txs(0,array(box,DIM));
 177     }
 178   }
 179   debug_gmx();
 180
 181   return bDone;
 182 }
 183
 184 static int recv_coords(t_commrec *cr,int *nviol,int *nit,int *k,
 185                        rvec x[],matrix box,bool bKeep,bool bDone)
 186 {
 187   int        nodeid,natom;
 188
 189   /* Check whether there is something from anyone */
 190   debug_gmx();
 191 #ifdef HAVE_MPI
 192   gmx_rx(MPI_ANY_SOURCE,record(nodeid));
 193   do {
 194     usleep(1000);
 195   } while (!mpiio_rx_probe(MPI_ANY_SOURCE));
 196 #endif
 197
 198   debug_gmx();
 199   if ((nodeid >= cr->nnodes) || (nodeid <= 0))
 200     fatal_error(0,"Reading data from nodeid %d",nodeid);
 201
 202   gmx_txs(nodeid,record(bDone));
 203
 204   if (!bDone) {
 205     gmx_rxs(nodeid,record(*nviol));
 206     gmx_rxs(nodeid,record(*nit));
 207     gmx_rxs(nodeid,record(*k));
 208     if (bKeep || (*nviol == 0)) {
 209       gmx_rxs(nodeid,record(natom));
 210       gmx_rxs(nodeid,array(x,natom));
 211       gmx_rxs(nodeid,array(box,DIM));
 212     }
 213   }
 214   return nodeid;
 215 }
 216
 217 void disco_slave(t_commrec *cr,FILE *fp)
 218 {
 219   t_correct *c;
 220   int       seed,nit,natom,nres,k,nviol;
 221   bool      bKeep,bDone;
 222   matrix    box;
 223   rvec      boxsize;
 224   rvec      *x,*xref,*xcenter;
 225
 226   debug_gmx();
 227   c = recv_init(fp,cr,&seed,&natom,&nres,&xref,&xcenter,&bKeep);
 228   c->nstprint = 0;
 229
 230   /* Make tags etc. */
 231   debug_gmx();
 232   init_corr2(c,natom);
 233
 234   debug_gmx();
 235   pr_corr(fp,c);
 236   fprintf(fp,"Random seed   = %d\n",seed);
 237
 238   snew(x,natom);
 239   bDone = FALSE;
 240   for(k=0; (!bDone); k++) {
 241     /* Generate random box*/
 242     debug_gmx();
 243     rand_box(c->bBox,box,boxsize,nres,c->bCubic,&seed);
 244
 245     /* Generate random coords */
 246     debug_gmx();
 247     rand_coords(natom,x,xref,c->weight,c->bCenter,xcenter,boxsize,&seed);
 248
 249     /* Now correct the random coords */
 250     debug_gmx();
 251     nviol = shake_coords(fp,FALSE,k,natom,xref,x,&seed,box,c,&nit);
 252
 253     debug_gmx();
 254     bDone = send_coords(cr,nviol,nit,k,natom,x,box,bKeep);
 255   }
 256 }
 257
 258 void disco_master(t_commrec *cr,FILE *fp,char *outfn,char *keepfn,t_correct *c,
 259                   bool bVerbose,t_atoms *atoms,
 260                   rvec xref[],rvec xcenter[],
 261                   int nstruct,int *seed,
 262                   bool bFit,int nfit,atom_id fit_ind[],
 263                   bool bPrintViol,char *violfn,rvec boxsize)
 264 {
 265   FILE    *gp;
 266   int     *nconvdist;
 267   int     i,k,kk,nconv,ntry,status,kstatus,natom,nres,nit;
 268   int     nvtest,nodeid,knodeid,nviol;
 269   double  tnit;
 270   rvec    *x,xcm;
 271   matrix  box,wrbox;
 272   atom_id *wr_ind;
 273   real    *w_rls;
 274   bool    bConverged;
 275
 276   natom = atoms->nr;
 277   nres  = atoms->nres;
 278
 279   /* Send out the word to my disciples */
 280   debug_gmx();
 281   send_init(fp,cr,c,seed,natom,nres,xref,xcenter,(keepfn != NULL));
 282
 283   /* Make tags etc. */
 284   debug_gmx();
 285   init_corr2(c,natom);
 286
 287   clear_mat(wrbox);
 288   wrbox[XX][XX] = wrbox[YY][YY] = wrbox[ZZ][ZZ] = nres;
 289   status = open_trx(outfn,"w");
 290   if (keepfn)
 291     kstatus = open_trx(keepfn,"w");
 292   else
 293     kstatus = -1;
 294   snew(x,natom);
 295   snew(wr_ind,natom);
 296   for(k=0; (k<natom); k++)
 297     wr_ind[k]=k;
 298   snew(w_rls,natom);
 299   for(k=0; (k<nfit); k++)
 300     w_rls[fit_ind[k]] = 1;
 301
 302   snew(nconvdist,c->maxnit+1);
 303   nconv = 0;
 304   ntry  = 0;
 305   tnit  = 0;
 306   debug_gmx();
 307   for(k=0; (k<nstruct); ) {
 308     nodeid        = recv_coords(cr,&nviol,&nit,&knodeid,x,box,
 309                                 (keepfn != NULL),FALSE);
 310     bConverged = (nviol == 0);
 311     tnit      += nit;
 312     ntry++;
 313
 314     if (bConverged)
 315       nconvdist[nit]++;
 316
 317     /*nvtest = quick_check(bVerbose ? fp : NULL,natom,x,box,c);*/
 318     if ((k % 20) == 0)
 319       fprintf(stderr,"%8s%6s%6s%8s%6s\n",
 320               "Struct","CPU","Str","nviol","niter");
 321     fprintf(stderr,"%8d%6d%6d%8d%6d\n",k,nodeid,knodeid,nviol,nit);
 322
 323     if (bConverged || keepfn) {
 324       center_in_box(natom,x,wrbox,x);
 325       if (bFit)
 326         do_fit(natom,w_rls,xref,x);
 327       write_trx(bConverged ? status : kstatus,
 328                 natom,wr_ind,atoms,k,(real) k,wrbox,x,NULL);
 329
 330       if (bConverged)
 331         nconv++;
 332
 333       k++;
 334     }
 335     if (bPrintViol && 0) {
 336       /* Print structure coloured by the violations */
 337       if (!atoms->pdbinfo)
 338         snew(atoms->pdbinfo,natom);
 339       for(kk=0; (kk<natom); kk++)
 340         atoms->pdbinfo[kk].bfac = (real) c->bViol[kk];
 341       gp=ffopen(violfn,"w");
 342       write_pdbfile(gp,"Structure coloured by violation",atoms,x,box,0,-1);
 343       ffclose(gp);
 344     }
 345   }
 346   for(k=1; (k<cr->nnodes); k++)
 347     nodeid = recv_coords(cr,&nviol,&nit,&knodeid,x,box,(keepfn != NULL),TRUE);
 348
 349   close_trx(status);
 350   if (keepfn)
 351     close_trx(kstatus);
 352   gp = xvgropen("conv_stat.xvg","Iterations per converged structure",
 353                 "nit","N");
 354   for(i=0; (i<c->maxnit); i++)
 355     fprintf(gp,"%10d  %10d\n",i,nconvdist[i]);
 356   ffclose(gp);
 357   sfree(x);
 358   sfree(w_rls);
 359   sfree(wr_ind);
 360   sfree(nconvdist);
 361
 362   pr_conv_stat(fp,ntry,nconv,tnit);
 363   pr_conv_stat(stderr,ntry,nconv,tnit);
 364 }
 365