src/tools/addconf.c

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