3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
52 #include "mtop_util.h"
53 #include "chargegroup.h"
55 static real box_margin;
57 static real max_dist(rvec *x, real *r, int start, int end)
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]));
70 static gmx_bool outside_box_minus_margin2(rvec x,matrix box)
72 return ( (x[XX]<2*box_margin) || (x[XX]>box[XX][XX]-2*box_margin) ||
73 (x[YY]<2*box_margin) || (x[YY]>box[YY][YY]-2*box_margin) ||
74 (x[ZZ]<2*box_margin) || (x[ZZ]>box[ZZ][ZZ]-2*box_margin) );
77 static gmx_bool outside_box_plus_margin(rvec x,matrix box)
79 return ( (x[XX]<-box_margin) || (x[XX]>box[XX][XX]+box_margin) ||
80 (x[YY]<-box_margin) || (x[YY]>box[YY][YY]+box_margin) ||
81 (x[ZZ]<-box_margin) || (x[ZZ]>box[ZZ][ZZ]+box_margin) );
84 static int mark_res(int at, gmx_bool *mark, int natoms, t_atom *atom,int *nmark)
88 resind = atom[at].resind;
89 while( (at > 0) && (resind==atom[at-1].resind) )
91 while( (at < natoms) && (resind==atom[at].resind) ) {
102 static real find_max_real(int n,real radius[])
111 rmax = max(rmax,radius[i]);
116 static void combine_atoms(t_atoms *ap,t_atoms *as,
117 rvec xp[],rvec *vp,rvec xs[],rvec *vs,
118 t_atoms **a_comb,rvec **x_comb,rvec **v_comb)
124 /* Total number of atoms */
125 natot = ap->nr+as->nr;
128 init_t_atoms(ac,natot,FALSE);
131 if (vp && vs) snew(vc,natot);
133 /* Fill the new structures */
134 for(i=j=0; (i<ap->nr); i++,j++) {
135 copy_rvec(xp[i],xc[j]);
136 if (vc) copy_rvec(vp[i],vc[j]);
137 memcpy(&(ac->atom[j]),&(ap->atom[i]),sizeof(ap->atom[i]));
138 ac->atom[j].type = 0;
141 for(i=0; (i<as->nr); i++,j++) {
142 copy_rvec(xs[i],xc[j]);
143 if (vc) copy_rvec(vs[i],vc[j]);
144 memcpy(&(ac->atom[j]),&(as->atom[i]),sizeof(as->atom[i]));
145 ac->atom[j].type = 0;
146 ac->atom[j].resind += res0;
149 ac->nres = ac->atom[j-1].resind+1;
150 /* Fill all elements to prevent uninitialized memory */
151 for(i=0; i<ac->nr; i++) {
156 ac->atom[i].type = 0;
157 ac->atom[i].typeB = 0;
158 ac->atom[i].ptype = eptAtom;
167 static t_forcerec *fr=NULL;
169 void do_nsgrid(FILE *fp,gmx_bool bVerbose,
170 matrix box,rvec x[],t_atoms *atoms,real rlong,
171 const output_env_t oenv)
190 /* Charge group index */
191 snew(cg_index,natoms);
192 for(i=0; (i<natoms); i++)
195 /* Topology needs charge groups and exclusions */
198 mtop->natoms = natoms;
199 /* Make one moltype that contains the whol system */
201 snew(mtop->moltype,mtop->nmoltype);
202 molt = &mtop->moltype[0];
203 molt->name = mtop->name;
204 molt->atoms = *atoms;
205 stupid_fill_block(&molt->cgs,mtop->natoms,FALSE);
206 stupid_fill_blocka(&molt->excls,natoms);
207 /* Make one molblock for the whole system */
209 snew(mtop->molblock,mtop->nmolblock);
210 mtop->molblock[0].type = 0;
211 mtop->molblock[0].nmol = 1;
212 mtop->molblock[0].natoms_mol = natoms;
213 /* Initialize a single energy group */
214 mtop->groups.grps[egcENER].nr = 1;
215 mtop->groups.ngrpnr[egcENER] = 0;
216 mtop->groups.grpnr[egcENER] = NULL;
218 ffp = &mtop->ffparams;
223 snew(ffp->functype,1);
224 snew(ffp->iparams,1);
225 ffp->iparams[0].lj.c6 = 1;
226 ffp->iparams[0].lj.c12 = 1;
228 /* inputrec structure */
230 ir->coulombtype = eelCUT;
231 ir->vdwtype = evdwCUT;
233 ir->ns_type = ensGRID;
234 snew(ir->opts.egp_flags,1);
236 top = gmx_mtop_generate_local_top(mtop,ir);
238 /* Some nasty shortcuts */
241 /* mdatoms structure */
244 md = init_mdatoms(fp,mtop,FALSE);
245 atoms2md(mtop,ir,0,NULL,0,mtop->natoms,md);
248 /* forcerec structure */
253 /* cr->nthreads = 1; */
255 /* ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
256 printf("Neighborsearching with a cut-off of %g\n",rlong);
257 init_forcerec(stdout,fr,ir,top,cr,md,box,FALSE,NULL,NULL,NULL,TRUE);*/
259 fr->hcg = top->cgs.nr;
262 /* Prepare for neighboursearching */
265 /* Init things dependent on parameters */
266 ir->rlistlong = ir->rlist = ir->rcoulomb = ir->rvdw = rlong;
267 /* create free energy data to avoid NULLs */
269 printf("Neighborsearching with a cut-off of %g\n",rlong);
270 init_forcerec(stdout,oenv,fr,NULL,ir,mtop,cr,box,FALSE,
271 NULL,NULL,NULL,NULL,NULL,TRUE,-1);
273 pr_forcerec(debug,fr,cr);
275 /* Calculate new stuff dependent on coords and box */
276 for(m=0; (m<DIM); m++)
277 box_size[m] = box[m][m];
278 calc_shifts(box,fr->shift_vec);
279 put_charge_groups_in_box(fp,0,cgs->nr,fr->ePBC,box,cgs,x,fr->cg_cm);
281 /* Do the actual neighboursearching */
284 init_neighbor_list(fp,fr,md->homenr);
285 search_neighbours(fp,fr,x,box,top,
286 &mtop->groups,cr,&nrnb,md,lambda,dvdl,NULL,TRUE,FALSE,FALSE);
289 dump_nblist(debug,cr,fr,0);
292 fprintf(stderr,"Successfully made neighbourlist\n");
295 gmx_bool bXor(gmx_bool b1,gmx_bool b2)
297 return (b1 && !b2) || (b2 && !b1);
300 void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, gmx_bool bSrenew,
301 int ePBC, matrix box, gmx_bool bInsert,
302 t_atoms *atoms_solvt,rvec *x_solvt,rvec *v_solvt,real *r_solvt,
303 gmx_bool bVerbose,real rshell,int max_sol, const output_env_t oenv)
307 real max_vdw,*r_prot,*r_all,n2,r2,ib1,ib2;
308 int natoms_prot,natoms_solvt;
309 int i,j,jj,m,j0,j1,jjj,jnres,jnr,inr,iprot,is1,is2;
310 int prev,resnr,nresadd,d,k,ncells,maxincell;
311 int dx0,dx1,dy0,dy1,dz0,dz1;
312 int ntest,nremove,nkeep;
313 rvec dx,xi,xj,xpp,*x_all,*v_all;
314 gmx_bool *remove,*keep;
317 natoms_prot = atoms->nr;
318 natoms_solvt = atoms_solvt->nr;
319 if (natoms_solvt <= 0) {
320 fprintf(stderr,"WARNING: Nothing to add\n");
324 if (ePBC == epbcSCREW)
325 gmx_fatal(FARGS,"Sorry, %s pbc is not yet supported",epbc_names[ePBC]);
328 fprintf(stderr,"Calculating Overlap...\n");
330 /* Set margin around box edges to largest solvent dimension.
331 * The maximum distance between atoms in a solvent molecule should
332 * be calculated. At the moment a fudge factor of 3 is used.
335 box_margin = 3*find_max_real(natoms_solvt,r_solvt);
336 max_vdw = max(3*find_max_real(natoms_prot,r_prot),box_margin);
337 fprintf(stderr,"box_margin = %g\n",box_margin);
339 snew(remove,natoms_solvt);
343 for(i=0; i<atoms_solvt->nr; i++)
344 if ( outside_box_plus_margin(x_solvt[i],box) )
345 i=mark_res(i,remove,atoms_solvt->nr,atoms_solvt->atom,&nremove);
346 fprintf(stderr,"Removed %d atoms that were outside the box\n",nremove);
349 /* Define grid stuff for genbox */
350 /* Largest VDW radius */
351 snew(r_all,natoms_prot+natoms_solvt);
352 for(i=j=0; i<natoms_prot; i++,j++)
354 for(i=0; i<natoms_solvt; i++,j++)
358 combine_atoms(atoms,atoms_solvt,*x,v?*v:NULL,x_solvt,v_solvt,
359 &atoms_all,&x_all,&v_all);
361 /* Do neighboursearching step */
362 do_nsgrid(stdout,bVerbose,box,x_all,atoms_all,max_vdw,oenv);
364 /* check solvent with solute */
365 nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
366 fprintf(stderr,"nri = %d, nrj = %d\n",nlist->nri,nlist->nrj);
367 for(bSolSol=0; (bSolSol<=(bInsert ? 0 : 1)); bSolSol++) {
369 fprintf(stderr,"Checking %s-Solvent overlap:",
370 bSolSol ? "Solvent" : "Protein");
371 for(i=0; (i<nlist->nri && nremove<natoms_solvt); i++) {
372 inr = nlist->iinr[i];
373 j0 = nlist->jindex[i];
374 j1 = nlist->jindex[i+1];
375 rvec_add(x_all[inr],fr->shift_vec[nlist->shift[i]],xi);
377 for(j=j0; (j<j1 && nremove<natoms_solvt); j++) {
378 jnr = nlist->jjnr[j];
379 copy_rvec(x_all[jnr],xj);
381 /* Check solvent-protein and solvent-solvent */
382 is1 = inr-natoms_prot;
383 is2 = jnr-natoms_prot;
385 /* Check if at least one of the atoms is a solvent that is not yet
386 * listed for removal, and if both are solvent, that they are not in the
390 bXor((is1 >= 0),(is2 >= 0)) && /* One atom is protein */
391 ((is1 < 0) || ((is1 >= 0) && !remove[is1])) &&
392 ((is2 < 0) || ((is2 >= 0) && !remove[is2]))) ||
395 (is1 >= 0) && (!remove[is1]) && /* is1 is solvent */
396 (is2 >= 0) && (!remove[is2]) && /* is2 is solvent */
397 (bInsert || /* when inserting also check inside the box */
398 (outside_box_minus_margin2(x_solvt[is1],box) && /* is1 on edge */
399 outside_box_minus_margin2(x_solvt[is2],box)) /* is2 on edge */
401 (atoms_solvt->atom[is1].resind != /* Not the same residue */
402 atoms_solvt->atom[is2].resind))) {
407 r2 = sqr(r_all[inr]+r_all[jnr]);
410 nremove = natoms_solvt;
411 for(k=0; k<nremove; k++) {
415 /* Need only remove one of the solvents... */
417 (void) mark_res(is2,remove,natoms_solvt,atoms_solvt->atom,
420 (void) mark_res(is1,remove,natoms_solvt,atoms_solvt->atom,
423 fprintf(stderr,"Neither atom is solvent%d %d\n",is1,is2);
429 fprintf(stderr," tested %d pairs, removed %d atoms.\n",ntest,nremove);
433 for(i=0; i<natoms_solvt; i++)
434 fprintf(debug,"remove[%5d] = %s\n",i,bool_names[remove[i]]);
436 /* Search again, now with another cut-off */
438 do_nsgrid(stdout,bVerbose,box,x_all,atoms_all,rshell,oenv);
439 nlist = &(fr->nblists[0].nlist_sr[eNL_VDW]);
440 fprintf(stderr,"nri = %d, nrj = %d\n",nlist->nri,nlist->nrj);
442 snew(keep,natoms_solvt);
443 for(i=0; i<nlist->nri; i++) {
444 inr = nlist->iinr[i];
445 j0 = nlist->jindex[i];
446 j1 = nlist->jindex[i+1];
448 for(j=j0; j<j1; j++) {
449 jnr = nlist->jjnr[j];
451 /* Check solvent-protein and solvent-solvent */
452 is1 = inr-natoms_prot;
453 is2 = jnr-natoms_prot;
455 /* Check if at least one of the atoms is a solvent that is not yet
456 * listed for removal, and if both are solvent, that they are not in the
460 mark_res(is1,keep,natoms_solvt,atoms_solvt->atom,&nkeep);
461 else if (is1<0 && is2>=0)
462 mark_res(is2,keep,natoms_solvt,atoms_solvt->atom,&nkeep);
465 fprintf(stderr,"Keeping %d solvent atoms after proximity check\n",
467 for (i=0; i<natoms_solvt; i++)
468 remove[i] = remove[i] || !keep[i];
471 /* count how many atoms and residues will be added and make space */
474 jnres = atoms_solvt->nres;
478 for (i=0; ((i<atoms_solvt->nr) &&
479 ((max_sol == 0) || (jnres < max_sol))); i++) {
483 (atoms_solvt->atom[i].resind != atoms_solvt->atom[i-1].resind))
489 fprintf(debug,"Will add %d atoms in %d residues\n",j,jnres);
491 /* Flag the remaing solvent atoms to be removed */
492 jjj = atoms_solvt->atom[i-1].resind;
493 for ( ; (i<atoms_solvt->nr); i++) {
494 if (atoms_solvt->atom[i].resind > jjj)
502 srenew(atoms->resinfo, atoms->nres+jnres);
503 srenew(atoms->atomname, atoms->nr+j);
504 srenew(atoms->atom, atoms->nr+j);
505 srenew(*x, atoms->nr+j);
506 if (v) srenew(*v, atoms->nr+j);
507 srenew(*r, atoms->nr+j);
510 /* add the selected atoms_solvt to atoms */
512 resnr = atoms->resinfo[atoms->atom[atoms->nr-1].resind].nr;
518 for (i=0; i<atoms_solvt->nr; i++) {
521 atoms_solvt->atom[i].resind != atoms_solvt->atom[prev].resind) {
525 atoms->resinfo[atoms->nres-1] =
526 atoms_solvt->resinfo[atoms_solvt->atom[i].resind];
527 atoms->resinfo[atoms->nres-1].nr = resnr;
528 /* calculate shift of the solvent molecule using the first atom */
529 copy_rvec(x_solvt[i],dx);
530 put_atoms_in_box(ePBC,box,1,&dx);
531 rvec_dec(dx,x_solvt[i]);
533 atoms->atom[atoms->nr] = atoms_solvt->atom[i];
534 atoms->atomname[atoms->nr] = atoms_solvt->atomname[i];
535 rvec_add(x_solvt[i],dx,(*x)[atoms->nr]);
536 if (v) copy_rvec(v_solvt[i],(*v)[atoms->nr]);
537 (*r)[atoms->nr] = r_solvt[i];
538 atoms->atom[atoms->nr].resind = atoms->nres-1;
544 srenew(atoms->resinfo, atoms->nres+nresadd);
547 fprintf(stderr,"Added %d molecules\n",nresadd);
550 done_atom(atoms_all);