src/tools/addconf.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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdlib.h>
  43 #include <string.h>
  44 #include "vec.h"
  45 #include "macros.h"
  46 #include "smalloc.h"
  47 #include "addconf.h"
  48 #include "force.h"
  49 #include "gstat.h"
  50 #include "names.h"
  51 #include "nsgrid.h"
  52 #include "mdatoms.h"
  53 #include "nrnb.h"
  54 #include "ns.h"
  55 #include "mtop_util.h"
  56 #include "chargegroup.h"
  57
  58 static real box_margin;
  59
  60 static real max_dist(rvec *x, real *r, int start, int end)
  61 {
  62   real maxd;
  63   int i,j;
  64
  65   maxd=0;
  66   for(i=start; i<end; i++)
  67     for(j=i+1; j<end; j++)
  68       maxd=max(maxd,sqrt(distance2(x[i],x[j]))+0.5*(r[i]+r[j]));
  69
  70   return 0.5*maxd;
  71 }
  72
  73 static gmx_bool outside_box_minus_margin2(rvec x,matrix box)
  74 {
  75   return ( (x[XX]<2*box_margin) || (x[XX]>box[XX][XX]-2*box_margin) ||
  76            (x[YY]<2*box_margin) || (x[YY]>box[YY][YY]-2*box_margin) ||
  77            (x[ZZ]<2*box_margin) || (x[ZZ]>box[ZZ][ZZ]-2*box_margin) );
  78 }
  79
  80 static gmx_bool outside_box_plus_margin(rvec x,matrix box)
  81 {
  82   return ( (x[XX]<-box_margin) || (x[XX]>box[XX][XX]+box_margin) ||
  83            (x[YY]<-box_margin) || (x[YY]>box[YY][YY]+box_margin) ||
  84            (x[ZZ]<-box_margin) || (x[ZZ]>box[ZZ][ZZ]+box_margin) );
  85 }
  86
  87 static int mark_res(int at, gmx_bool *mark, int natoms, t_atom *atom,int *nmark)
  88 {
  89   int resind;
  90
  91   resind = atom[at].resind;
  92   while( (at > 0) && (resind==atom[at-1].resind) )
  93     at--;
  94   while( (at < natoms) && (resind==atom[at].resind) ) {
  95     if (!mark[at]) {
  96       mark[at]=TRUE;
  97       (*nmark)++;
  98     }
  99     at++;
 100   }
 101
 102   return at;
 103 }
 104
 105 static real find_max_real(int n,real radius[])
 106 {
 107   int  i;
 108   real rmax;
 109
 110   rmax = 0;
 111   if (n > 0) {
 112     rmax = radius[0];
 113     for(i=1; (i<n); i++)
 114       rmax = max(rmax,radius[i]);
 115   }
 116   return rmax;
 117 }
 118
 119 static void combine_atoms(t_atoms *ap,t_atoms *as,
 120                           rvec xp[],rvec *vp,rvec xs[],rvec *vs,
 121                           t_atoms **a_comb,rvec **x_comb,rvec **v_comb)
 122 {
 123   t_atoms *ac;
 124   rvec    *xc,*vc=NULL;
 125   int     i,j,natot,res0;
 126
 127   /* Total number of atoms */
 128   natot = ap->nr+as->nr;
 129
 130   snew(ac,1);
 131   init_t_atoms(ac,natot,FALSE);
 132
 133   snew(xc,natot);
 134   if (vp && vs) snew(vc,natot);
 135
 136   /* Fill the new structures */
 137   for(i=j=0; (i<ap->nr); i++,j++) {
 138     copy_rvec(xp[i],xc[j]);
 139     if (vc) copy_rvec(vp[i],vc[j]);
 140     memcpy(&(ac->atom[j]),&(ap->atom[i]),sizeof(ap->atom[i]));
 141     ac->atom[j].type = 0;
 142   }
 143   res0 = ap->nres;
 144   for(i=0; (i<as->nr); i++,j++) {
 145     copy_rvec(xs[i],xc[j]);
 146     if (vc) copy_rvec(vs[i],vc[j]);
 147     memcpy(&(ac->atom[j]),&(as->atom[i]),sizeof(as->atom[i]));
 148     ac->atom[j].type    = 0;
 149     ac->atom[j].resind += res0;
 150   }
 151   ac->nr   = j;
 152   ac->nres = ac->atom[j-1].resind+1;
 153   /* Fill all elements to prevent uninitialized memory */
 154   for(i=0; i<ac->nr; i++) {
 155     ac->atom[i].m     = 1;
 156     ac->atom[i].q     = 0;
 157     ac->atom[i].mB    = 1;
 158     ac->atom[i].qB    = 0;
 159     ac->atom[i].type  = 0;
 160     ac->atom[i].typeB = 0;
 161     ac->atom[i].ptype = eptAtom;
 162   }
 163
 164   /* Return values */
 165   *a_comb = ac;
 166   *x_comb = xc;
 167   *v_comb = vc;
 168 }
 169
 170 static t_forcerec *fr=NULL;
 171
 172 void do_nsgrid(FILE *fp,gmx_bool bVerbose,
 173                matrix box,rvec x[],t_atoms *atoms,real rlong,
 174                const output_env_t oenv)
 175 {
 176   gmx_mtop_t *mtop;
 177   gmx_localtop_t *top;
 178   t_mdatoms  *md;
 179   t_block    *cgs;
 180   t_inputrec *ir;
 181   t_nrnb     nrnb;
 182   t_commrec  *cr;
 183   int        *cg_index;
 184   gmx_moltype_t *molt;
 185   gmx_ffparams_t *ffp;
 186   ivec       *nFreeze;
 187   int        i,m,natoms;
 188   rvec       box_size;
 189   real       *lambda,*dvdl;
 190
 191   natoms = atoms->nr;
 192
 193   /* Charge group index */
 194   snew(cg_index,natoms);
 195   for(i=0; (i<natoms); i++)
 196     cg_index[i]=i;
 197
 198   /* Topology needs charge groups and exclusions */
 199   snew(mtop,1);
 200   init_mtop(mtop);
 201   mtop->natoms = natoms;
 202   /* Make one moltype that contains the whol system */
 203   mtop->nmoltype = 1;
 204   snew(mtop->moltype,mtop->nmoltype);
 205   molt = &mtop->moltype[0];
 206   molt->name  = mtop->name;
 207   molt->atoms = *atoms;
 208   stupid_fill_block(&molt->cgs,mtop->natoms,FALSE);
 209   stupid_fill_blocka(&molt->excls,natoms);
 210   /* Make one molblock for the whole system */
 211   mtop->nmolblock = 1;
 212   snew(mtop->molblock,mtop->nmolblock);
 213   mtop->molblock[0].type = 0;
 214   mtop->molblock[0].nmol = 1;
 215   mtop->molblock[0].natoms_mol = natoms;
 216   /* Initialize a single energy group */
 217   mtop->groups.grps[egcENER].nr = 1;
 218   mtop->groups.ngrpnr[egcENER]  = 0;
 219   mtop->groups.grpnr[egcENER]   = NULL;
 220
 221   ffp = &mtop->ffparams;
 222
 223   ffp->ntypes = 1;
 224   ffp->atnr   = 1;
 225   ffp->reppow = 12;
 226   snew(ffp->functype,1);
 227   snew(ffp->iparams,1);
 228   ffp->iparams[0].lj.c6  = 1;
 229   ffp->iparams[0].lj.c12 = 1;
 230
 231   /* inputrec structure */
 232   snew(ir,1);
 233   ir->coulombtype = eelCUT;
 234   ir->vdwtype     = evdwCUT;
 235   ir->ndelta      = 2;
 236   ir->ns_type     = ensGRID;
 237   snew(ir->opts.egp_flags,1);
 238
 239   top = gmx_mtop_generate_local_top(mtop,ir);
 240
 241   /* Some nasty shortcuts */
 242   cgs  = &(top->cgs);
 243
 244     /* mdatoms structure */
 245   snew(nFreeze,2);
 246   snew(md,1);
 247   md = init_mdatoms(fp,mtop,FALSE);
 248   atoms2md(mtop,ir,0,NULL,0,mtop->natoms,md);
 249   sfree(nFreeze);
 250
 251   /* forcerec structure */
 252   if (fr == NULL)
 253     fr = mk_forcerec();
 254   snew(cr,1);
 255   cr->nnodes   = 1;
 256   /* cr->nthreads = 1; */
 257
 258   /*    ir->rlist       = ir->rcoulomb = ir->rvdw = rlong;
 259         printf("Neighborsearching with a cut-off of %g\n",rlong);
 260         init_forcerec(stdout,fr,ir,top,cr,md,box,FALSE,NULL,NULL,NULL,TRUE);*/
 261   fr->cg0 = 0;
 262   fr->hcg = top->cgs.nr;
 263   fr->nWatMol = 0;
 264
 265   /* Prepare for neighboursearching */
 266   init_nrnb(&nrnb);
 267
 268   /* Init things dependent on parameters */
 269   ir->rlistlong = ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
 270   /* create free energy data to avoid NULLs */
 271   snew(ir->fepvals,1);
 272   printf("Neighborsearching with a cut-off of %g\n",rlong);
 273   init_forcerec(stdout,oenv,fr,NULL,ir,mtop,cr,box,FALSE,
 274                 NULL,NULL,NULL,NULL,NULL,TRUE,-1);
 275   if (debug)
 276     pr_forcerec(debug,fr,cr);
 277
 278   /* Calculate new stuff dependent on coords and box */
 279   for(m=0; (m<DIM); m++)
 280     box_size[m] = box[m][m];
 281   calc_shifts(box,fr->shift_vec);
 282   put_charge_groups_in_box(fp,0,cgs->nr,fr->ePBC,box,cgs,x,fr->cg_cm);
 283
 284   /* Do the actual neighboursearching */
 285   snew(lambda,efptNR);
 286   snew(dvdl,efptNR);
 287   init_neighbor_list(fp,fr,md->homenr);
 288   search_neighbours(fp,fr,x,box,top,
 289                     &mtop->groups,cr,&nrnb,md,lambda,dvdl,NULL,TRUE,FALSE,FALSE);
 290
 291   if (debug)
 292     dump_nblist(debug,cr,fr,0);
 293
 294   if (bVerbose)
 295     fprintf(stderr,"Successfully made neighbourlist\n");
 296 }
 297
 298 gmx_bool bXor(gmx_bool b1,gmx_bool b2)
 299 {
 300   return (b1 && !b2) || (b2 && !b1);
 301 }
 302
 303 void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, gmx_bool bSrenew,
 304               int ePBC, matrix box, gmx_bool bInsert,
 305               t_atoms *atoms_solvt,rvec *x_solvt,rvec *v_solvt,real *r_solvt,
 306               gmx_bool bVerbose,real rshell,int max_sol, const output_env_t oenv)
 307 {
 308   t_nblist   *nlist;
 309   t_atoms    *atoms_all;
 310   real       max_vdw,*r_prot,*r_all,n2,r2,ib1,ib2;
 311   int        natoms_prot,natoms_solvt;
 312   int        i,j,jj,m,j0,j1,jjj,jnres,jnr,inr,iprot,is1,is2;
 313   int        prev,resnr,nresadd,d,k,ncells,maxincell;
 314   int        dx0,dx1,dy0,dy1,dz0,dz1;
 315   int        ntest,nremove,nkeep;
 316   rvec       dx,xi,xj,xpp,*x_all,*v_all;
 317   gmx_bool       *remove,*keep;
 318   int        bSolSol;
 319
 320   natoms_prot  = atoms->nr;
 321   natoms_solvt = atoms_solvt->nr;
 322   if (natoms_solvt <= 0) {
 323     fprintf(stderr,"WARNING: Nothing to add\n");
 324     return;
 325   }
 326
 327   if (ePBC == epbcSCREW)
 328     gmx_fatal(FARGS,"Sorry, %s pbc is not yet supported",epbc_names[ePBC]);
 329
 330   if (bVerbose)
 331     fprintf(stderr,"Calculating Overlap...\n");
 332
 333   /* Set margin around box edges to largest solvent dimension.
 334    * The maximum distance between atoms in a solvent molecule should
 335    * be calculated. At the moment a fudge factor of 3 is used.
 336    */
 337   r_prot     = *r;
 338   box_margin = 3*find_max_real(natoms_solvt,r_solvt);
 339   max_vdw    = max(3*find_max_real(natoms_prot,r_prot),box_margin);
 340   fprintf(stderr,"box_margin = %g\n",box_margin);
 341
 342   snew(remove,natoms_solvt);
 343
 344   nremove = 0;
 345   if (!bInsert) {
 346     for(i=0; i<atoms_solvt->nr; i++)
 347       if ( outside_box_plus_margin(x_solvt[i],box) )
 348         i=mark_res(i,remove,atoms_solvt->nr,atoms_solvt->atom,&nremove);
 349     fprintf(stderr,"Removed %d atoms that were outside the box\n",nremove);
 350   }
 351
 352   /* Define grid stuff for genbox */
 353   /* Largest VDW radius */
 354   snew(r_all,natoms_prot+natoms_solvt);
 355   for(i=j=0; i<natoms_prot; i++,j++)
 356     r_all[j]=r_prot[i];
 357   for(i=0; i<natoms_solvt; i++,j++)
 358     r_all[j]=r_solvt[i];
 359
 360   /* Combine arrays */
 361   combine_atoms(atoms,atoms_solvt,*x,v?*v:NULL,x_solvt,v_solvt,
 362                 &atoms_all,&x_all,&v_all);
 363
 364   /* Do neighboursearching step */
 365   do_nsgrid(stdout,bVerbose,box,x_all,atoms_all,max_vdw,oenv);
 366
 367   /* check solvent with solute */
 368   nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
 369   fprintf(stderr,"nri = %d, nrj = %d\n",nlist->nri,nlist->nrj);
 370   for(bSolSol=0; (bSolSol<=(bInsert ? 0 : 1)); bSolSol++) {
 371     ntest = nremove = 0;
 372     fprintf(stderr,"Checking %s-Solvent overlap:",
 373             bSolSol ? "Solvent" : "Protein");
 374     for(i=0; (i<nlist->nri && nremove<natoms_solvt); i++) {
 375       inr = nlist->iinr[i];
 376       j0  = nlist->jindex[i];
 377       j1  = nlist->jindex[i+1];
 378       rvec_add(x_all[inr],fr->shift_vec[nlist->shift[i]],xi);
 379
 380       for(j=j0; (j<j1 && nremove<natoms_solvt); j++) {
 381         jnr = nlist->jjnr[j];
 382         copy_rvec(x_all[jnr],xj);
 383
 384         /* Check solvent-protein and solvent-solvent */
 385         is1 = inr-natoms_prot;
 386         is2 = jnr-natoms_prot;
 387
 388         /* Check if at least one of the atoms is a solvent that is not yet
 389          * listed for removal, and if both are solvent, that they are not in the
 390          * same residue.
 391          */
 392         if ((!bSolSol &&
 393              bXor((is1 >= 0),(is2 >= 0)) &&  /* One atom is protein */
 394              ((is1 < 0) || ((is1 >= 0) && !remove[is1])) &&
 395              ((is2 < 0) || ((is2 >= 0) && !remove[is2]))) ||
 396
 397             (bSolSol  &&
 398              (is1 >= 0) && (!remove[is1]) &&   /* is1 is solvent */
 399              (is2 >= 0) && (!remove[is2]) &&   /* is2 is solvent */
 400              (bInsert || /* when inserting also check inside the box */
 401               (outside_box_minus_margin2(x_solvt[is1],box) && /* is1 on edge */
 402                outside_box_minus_margin2(x_solvt[is2],box))   /* is2 on edge */
 403               ) &&
 404              (atoms_solvt->atom[is1].resind !=  /* Not the same residue */
 405               atoms_solvt->atom[is2].resind))) {
 406
 407           ntest++;
 408           rvec_sub(xi,xj,dx);
 409           n2 = norm2(dx);
 410           r2 = sqr(r_all[inr]+r_all[jnr]);
 411           if (n2 < r2) {
 412             if (bInsert) {
 413               nremove = natoms_solvt;
 414               for(k=0; k<nremove; k++) {
 415                 remove[k] = TRUE;
 416               }
 417             }
 418             /* Need only remove one of the solvents... */
 419             if (is2 >= 0)
 420               (void) mark_res(is2,remove,natoms_solvt,atoms_solvt->atom,
 421                               &nremove);
 422             else if (is1 >= 0)
 423               (void) mark_res(is1,remove,natoms_solvt,atoms_solvt->atom,
 424                               &nremove);
 425             else
 426               fprintf(stderr,"Neither atom is solvent%d %d\n",is1,is2);
 427           }
 428         }
 429       }
 430     }
 431     if (!bInsert) {
 432       fprintf(stderr," tested %d pairs, removed %d atoms.\n",ntest,nremove);
 433     }
 434   }
 435   if (debug)
 436     for(i=0; i<natoms_solvt; i++)
 437       fprintf(debug,"remove[%5d] = %s\n",i,bool_names[remove[i]]);
 438
 439   /* Search again, now with another cut-off */
 440   if (rshell > 0) {
 441     do_nsgrid(stdout,bVerbose,box,x_all,atoms_all,rshell,oenv);
 442     nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
 443     fprintf(stderr,"nri = %d, nrj = %d\n",nlist->nri,nlist->nrj);
 444     nkeep = 0;
 445     snew(keep,natoms_solvt);
 446     for(i=0; i<nlist->nri; i++) {
 447       inr = nlist->iinr[i];
 448       j0  = nlist->jindex[i];
 449       j1  = nlist->jindex[i+1];
 450
 451       for(j=j0; j<j1; j++) {
 452         jnr = nlist->jjnr[j];
 453
 454         /* Check solvent-protein and solvent-solvent */
 455         is1 = inr-natoms_prot;
 456         is2 = jnr-natoms_prot;
 457
 458         /* Check if at least one of the atoms is a solvent that is not yet
 459          * listed for removal, and if both are solvent, that they are not in the
 460          * same residue.
 461          */
 462         if (is1>=0 && is2<0)
 463           mark_res(is1,keep,natoms_solvt,atoms_solvt->atom,&nkeep);
 464         else if (is1<0 && is2>=0)
 465           mark_res(is2,keep,natoms_solvt,atoms_solvt->atom,&nkeep);
 466       }
 467     }
 468     fprintf(stderr,"Keeping %d solvent atoms after proximity check\n",
 469             nkeep);
 470     for (i=0; i<natoms_solvt; i++)
 471       remove[i] = remove[i] || !keep[i];
 472     sfree(keep);
 473   }
 474   /* count how many atoms and residues will be added and make space */
 475   if (bInsert) {
 476     j     = atoms_solvt->nr;
 477     jnres = atoms_solvt->nres;
 478   } else {
 479     j     = 0;
 480     jnres = 0;
 481     for (i=0; ((i<atoms_solvt->nr) &&
 482                ((max_sol == 0) || (jnres < max_sol))); i++) {
 483       if (!remove[i]) {
 484         j++;
 485         if ((i == 0) ||
 486             (atoms_solvt->atom[i].resind != atoms_solvt->atom[i-1].resind))
 487           jnres++;
 488       }
 489     }
 490   }
 491   if (debug)
 492     fprintf(debug,"Will add %d atoms in %d residues\n",j,jnres);
 493   if (!bInsert) {
 494     /* Flag the remaing solvent atoms to be removed */
 495     jjj = atoms_solvt->atom[i-1].resind;
 496     for ( ; (i<atoms_solvt->nr); i++) {
 497       if (atoms_solvt->atom[i].resind > jjj)
 498         remove[i] = TRUE;
 499       else
 500         j++;
 501     }
 502   }
 503
 504   if (bSrenew) {
 505     srenew(atoms->resinfo,  atoms->nres+jnres);
 506     srenew(atoms->atomname, atoms->nr+j);
 507     srenew(atoms->atom,     atoms->nr+j);
 508     srenew(*x,              atoms->nr+j);
 509     if (v) srenew(*v,       atoms->nr+j);
 510     srenew(*r,              atoms->nr+j);
 511   }
 512
 513   /* add the selected atoms_solvt to atoms */
 514   if (atoms->nr > 0) {
 515     resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
 516   } else {
 517     resnr = 0;
 518   }
 519   prev = -1;
 520   nresadd = 0;
 521   for (i=0; i<atoms_solvt->nr; i++) {
 522     if (!remove[i]) {
 523       if (prev == -1 ||
 524           atoms_solvt->atom[i].resind != atoms_solvt->atom[prev].resind) {
 525         nresadd ++;
 526         atoms->nres++;
 527         resnr++;
 528         atoms->resinfo[atoms->nres-1] =
 529           atoms_solvt->resinfo[atoms_solvt->atom[i].resind];
 530         atoms->resinfo[atoms->nres-1].nr = resnr;
 531         /* calculate shift of the solvent molecule using the first atom */
 532         copy_rvec(x_solvt[i],dx);
 533         put_atoms_in_box(ePBC,box,1,&dx);
 534         rvec_dec(dx,x_solvt[i]);
 535       }
 536       atoms->atom[atoms->nr] = atoms_solvt->atom[i];
 537       atoms->atomname[atoms->nr] = atoms_solvt->atomname[i];
 538       rvec_add(x_solvt[i],dx,(*x)[atoms->nr]);
 539       if (v) copy_rvec(v_solvt[i],(*v)[atoms->nr]);
 540       (*r)[atoms->nr]   = r_solvt[i];
 541       atoms->atom[atoms->nr].resind = atoms->nres-1;
 542       atoms->nr++;
 543       prev=i;
 544     }
 545   }
 546   if (bSrenew)
 547     srenew(atoms->resinfo,  atoms->nres+nresadd);
 548
 549   if (bVerbose)
 550     fprintf(stderr,"Added %d molecules\n",nresadd);
 551
 552   sfree(remove);
 553   done_atom(atoms_all);
 554   sfree(x_all);
 555   sfree(v_all);
 556 }