src/tools/addconf.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_addconf_c = "$Id$";
  30
  31 #include "vec.h"
  32 #include "macros.h"
  33 #include "smalloc.h"
  34 #include "addconf.h"
  35 #include "gstat.h"
  36 #include "princ.h"
  37 #include "rdgroup.h"
  38 #include "txtdump.h"
  39 #include "pbc.h"
  40
  41 static real box_margin;
  42
  43 real max_dist(rvec *x, real *r, int start, int end)
  44 {
  45   real maxd;
  46   int i,j;
  47
  48   maxd=0;
  49   for(i=start; i<end; i++)
  50     for(j=i+1; j<end; j++)
  51       maxd=max(maxd,sqrt(distance2(x[i],x[j]))+0.5*(r[i]+r[j]));
  52
  53   return 0.5*maxd;
  54 }
  55
  56 void set_margin(t_atoms *atoms, rvec *x, real *r)
  57 {
  58   int i,d,start;
  59 /*   char *resname; */
  60
  61   box_margin=0;
  62
  63   start=0;
  64   for(i=0; i < atoms->nr; i++) {
  65     if ( (i+1 == atoms->nr) ||
  66          (atoms->atom[i+1].resnr != atoms->atom[i].resnr) ) {
  67       d=max_dist(x,r,start,i+1);
  68       if (debug && d>box_margin)
  69         fprintf(debug,"getting margin from %s: %g\n",
  70                 *(atoms->resname[atoms->atom[i].resnr]),box_margin);
  71       box_margin=max(box_margin,d);
  72       start=i+1;
  73     }
  74   }
  75 }
  76
  77
  78 bool in_box_plus_margin(rvec x,matrix box)
  79 {
  80   return ( x[XX]>=-box_margin && x[XX]<=box[XX][XX]+box_margin &&
  81            x[YY]>=-box_margin && x[YY]<=box[YY][YY]+box_margin &&
  82            x[ZZ]>=-box_margin && x[ZZ]<=box[ZZ][ZZ]+box_margin );
  83 }
  84
  85 bool outside_box_minus_margin(rvec x,matrix box)
  86 {
  87   return ( x[XX]<box_margin || x[XX]>box[XX][XX]-box_margin ||
  88            x[YY]<box_margin || x[YY]>box[YY][YY]-box_margin ||
  89            x[ZZ]<box_margin || x[ZZ]>box[ZZ][ZZ]-box_margin );
  90 }
  91
  92 int mark_remove_res(int at, bool *remove, int natoms, t_atom *atom)
  93 {
  94   int resnr;
  95
  96   resnr = atom[at].resnr;
  97   while( (at > 0) && (resnr==atom[at-1].resnr) )
  98     at--;
  99   while( (at < natoms) && (resnr==atom[at].resnr) ) {
 100     remove[at]=TRUE;
 101     at++;
 102   }
 103
 104   return at;
 105 }
 106
 107 void add_conf(t_atoms *atoms, rvec **x, real **r,  bool bSrenew,
 108               int NTB, matrix box,
 109               t_atoms *atoms_solvt, rvec *x_solvt, real *r_solvt,
 110               bool bVerbose)
 111 {
 112   int  i,j,m,prev,resnr,nresadd,d,k;
 113   int  *atom_flag;
 114   rvec dx;
 115   bool *remove;
 116
 117   if (atoms_solvt->nr <= 0) {
 118     fprintf(stderr,"WARNING: Nothing to add\n");
 119     return;
 120   }
 121
 122   if (bVerbose)
 123     fprintf(stderr,"Calculating Overlap...\n");
 124
 125   snew(remove,atoms_solvt->nr);
 126   init_pbc(box,FALSE);
 127
 128   /* set margin around box edges to largest solvent dimension */
 129   set_margin(atoms_solvt,x_solvt,r_solvt);
 130
 131   /* remove atoms that are far outside the box */
 132   for(i=0; i<atoms_solvt->nr ;i++)
 133     if ( !in_box_plus_margin(x_solvt[i],box) )
 134       i=mark_remove_res(i,remove,atoms_solvt->nr,atoms_solvt->atom);
 135
 136   /* check solvent with solute */
 137   for(i=0; i<atoms->nr; i++) {
 138     if ( bVerbose && ( (i<10) || ((i+1)%10) || (i>atoms->nr-10) ) )
 139       fprintf(stderr,"\r%d out of %d atoms checked",i+1,atoms->nr);
 140     for(j=0; j<atoms_solvt->nr; j++)
 141       /* only check solvent that hasn't been marked for removal already */
 142       if (!remove[j]) {
 143         pbc_dx((*x)[i],x_solvt[j],dx);
 144         if ( norm2(dx) < sqr((*r)[i] + r_solvt[j]) )
 145           j=mark_remove_res(j,remove,atoms_solvt->nr,atoms_solvt->atom);
 146       }
 147   }
 148   if (bVerbose)
 149     fprintf(stderr,"\n");
 150
 151   /* check solvent with itself */
 152   for(i=0; i<atoms_solvt->nr ;i++) {
 153     if ( bVerbose && ( (i<10) || !((i+1)%10) || (i>atoms_solvt->nr-10) ) )
 154       fprintf(stderr,"\rchecking atom %d out of %d",i+1,atoms_solvt->nr);
 155
 156     /* check only atoms that haven't been marked for removal already and
 157        are close to the box edges */
 158     if ( !remove[i] && outside_box_minus_margin(x_solvt[i],box) )
 159       /* check with other border atoms,
 160          only if not already removed and two different molecules (resnr) */
 161       for(j=i+1; (j<atoms_solvt->nr) && !remove[i]; j++)
 162         if ( !remove[j] &&
 163              atoms_solvt->atom[i].resnr != atoms_solvt->atom[j].resnr &&
 164              outside_box_minus_margin(x_solvt[j],box) ) {
 165           pbc_dx(x_solvt[i],x_solvt[j],dx);
 166           if ( norm2(dx) < sqr(r_solvt[i] + r_solvt[j]) )
 167             j=mark_remove_res(j,remove,atoms_solvt->nr,atoms_solvt->atom);
 168         }
 169   }
 170   if (bVerbose)
 171     fprintf(stderr,"\n");
 172
 173   /* count how many atoms will be added and make space */
 174   j=0;
 175   for (i=0; (i<atoms_solvt->nr); i++)
 176     if (!remove[i])
 177       j++;
 178   if (bSrenew) {
 179     srenew(atoms->resname,  atoms->nres+atoms_solvt->nres);
 180     srenew(atoms->atomname, atoms->nr+j);
 181     srenew(atoms->atom,     atoms->nr+j);
 182     srenew(*x,              atoms->nr+j);
 183     srenew(*r,              atoms->nr+j);
 184   }
 185
 186   /* add the selected atoms_solvt to atoms */
 187   prev=NOTSET;
 188   nresadd=0;
 189   for (i=0; (i<atoms_solvt->nr); i++)
 190     if (!remove[i]) {
 191       if ( (prev==NOTSET) ||
 192            (atoms_solvt->atom[i].resnr != atoms_solvt->atom[prev].resnr) ) {
 193         nresadd ++;
 194         atoms->nres++;
 195       }
 196       atoms->nr++;
 197       atoms->atomname[atoms->nr-1] = atoms_solvt->atomname[i];
 198       copy_rvec(x_solvt[i],(*x)[atoms->nr-1]);
 199       (*r)[atoms->nr-1]   = r_solvt[i];
 200       atoms->atom[atoms->nr-1].resnr = atoms->nres-1;
 201       atoms->resname[atoms->nres-1] =
 202         atoms_solvt->resname[atoms_solvt->atom[i].resnr];
 203       prev=i;
 204     }
 205   if (bSrenew)
 206     srenew(atoms->resname,  atoms->nres+nresadd);
 207
 208   if (bVerbose)
 209     fprintf(stderr,"Added %d molecules\n",nresadd);
 210
 211   sfree(remove);
 212 }
 213
 214 void orient_mol(t_atoms *atoms,char *indexnm,rvec x[], rvec *v)
 215 {
 216   int     isize;
 217   atom_id *index,*simp;
 218   char    *grpnames;
 219   real    totmass;
 220   int     i,m;
 221   rvec    xcm,prcomp;
 222   matrix  trans;
 223
 224   /* Make an index for principal component analysis */
 225   fprintf(stderr,"\nSelect group for orientation of molecule:\n");
 226   get_index(atoms,indexnm,1,&isize,&index,&grpnames);
 227   snew(simp,atoms->nr);
 228   for(i=0; (i<atoms->nr); i++) {
 229     simp[i]=i;
 230     atoms->atom[i].m=1;
 231   }
 232   totmass = sub_xcm(x,atoms->nr,simp,atoms->atom,xcm,FALSE);
 233   principal_comp(isize,index,atoms->atom,x,trans,prcomp);
 234
 235   /* Check whether this trans matrix mirrors the molecule */
 236   if (det(trans) < 0) {
 237     if (debug)
 238       fprintf(stderr,"Mirroring rotation matrix in Z direction\n");
 239     for(m=0; (m<DIM); m++)
 240       trans[ZZ][m] = -trans[ZZ][m];
 241   }
 242   rotate_atoms(atoms->nr,simp,x,trans);
 243
 244   if (debug) {
 245     pr_rvecs(stderr,0,"Rot Matrix",trans,DIM);
 246     fprintf(stderr,"Det(trans) = %g\n",det(trans));
 247
 248     /* print principal component data */
 249     fprintf(stderr,"Norm of principal axes before rotation: "
 250             "(%.3f, %.3f, %.3f)\n",prcomp[XX],prcomp[YY],prcomp[ZZ]);
 251     fprintf(stderr,"Totmass = %g\n",totmass);
 252     principal_comp(isize,index,atoms->atom,x,trans,prcomp);
 253     rotate_atoms(atoms->nr,simp,x,trans);
 254     if (v)
 255       rotate_atoms(atoms->nr,simp,v,trans);
 256     pr_rvecs(stderr,0,"Rot Matrix",trans,DIM);
 257   }
 258   sfree(simp);
 259   sfree(index);
 260   sfree(grpnames);
 261 }
 262