2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
55 #include "mtop_util.h"
56 #include "chargegroup.h"
58 static real box_margin
;
60 static real
max_dist(rvec
*x
, real
*r
, int start
, int end
)
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
]));
73 static gmx_bool
outside_box_minus_margin2(rvec x
,matrix box
)
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
) );
80 static gmx_bool
outside_box_plus_margin(rvec x
,matrix box
)
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
) );
87 static int mark_res(int at
, gmx_bool
*mark
, int natoms
, t_atom
*atom
,int *nmark
)
91 resind
= atom
[at
].resind
;
92 while( (at
> 0) && (resind
==atom
[at
-1].resind
) )
94 while( (at
< natoms
) && (resind
==atom
[at
].resind
) ) {
105 static real
find_max_real(int n
,real radius
[])
114 rmax
= max(rmax
,radius
[i
]);
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
)
127 /* Total number of atoms */
128 natot
= ap
->nr
+as
->nr
;
131 init_t_atoms(ac
,natot
,FALSE
);
134 if (vp
&& vs
) snew(vc
,natot
);
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;
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
;
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
++) {
159 ac
->atom
[i
].type
= 0;
160 ac
->atom
[i
].typeB
= 0;
161 ac
->atom
[i
].ptype
= eptAtom
;
170 static t_forcerec
*fr
=NULL
;
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
)
193 /* Charge group index */
194 snew(cg_index
,natoms
);
195 for(i
=0; (i
<natoms
); i
++)
198 /* Topology needs charge groups and exclusions */
201 mtop
->natoms
= natoms
;
202 /* Make one moltype that contains the whol system */
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 */
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
;
221 ffp
= &mtop
->ffparams
;
226 snew(ffp
->functype
,1);
227 snew(ffp
->iparams
,1);
228 ffp
->iparams
[0].lj
.c6
= 1;
229 ffp
->iparams
[0].lj
.c12
= 1;
231 /* inputrec structure */
233 ir
->coulombtype
= eelCUT
;
234 ir
->vdwtype
= evdwCUT
;
236 ir
->ns_type
= ensGRID
;
237 snew(ir
->opts
.egp_flags
,1);
239 top
= gmx_mtop_generate_local_top(mtop
,ir
);
241 /* Some nasty shortcuts */
244 /* mdatoms structure */
247 md
= init_mdatoms(fp
,mtop
,FALSE
);
248 atoms2md(mtop
,ir
,0,NULL
,0,mtop
->natoms
,md
);
251 /* forcerec structure */
256 /* cr->nthreads = 1; */
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);*/
262 fr
->hcg
= top
->cgs
.nr
;
265 /* Prepare for neighboursearching */
268 /* Init things dependent on parameters */
269 ir
->rlistlong
= ir
->rlist
= ir
->rcoulomb
= ir
->rvdw
= rlong
;
270 /* create free energy data to avoid NULLs */
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);
276 pr_forcerec(debug
,fr
,cr
);
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
);
284 /* Do the actual neighboursearching */
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
);
292 dump_nblist(debug
,cr
,fr
,0);
295 fprintf(stderr
,"Successfully made neighbourlist\n");
298 gmx_bool
bXor(gmx_bool b1
,gmx_bool b2
)
300 return (b1
&& !b2
) || (b2
&& !b1
);
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
)
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
;
320 natoms_prot
= atoms
->nr
;
321 natoms_solvt
= atoms_solvt
->nr
;
322 if (natoms_solvt
<= 0) {
323 fprintf(stderr
,"WARNING: Nothing to add\n");
327 if (ePBC
== epbcSCREW
)
328 gmx_fatal(FARGS
,"Sorry, %s pbc is not yet supported",epbc_names
[ePBC
]);
331 fprintf(stderr
,"Calculating Overlap...\n");
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.
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
);
342 snew(remove
,natoms_solvt
);
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
);
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
++)
357 for(i
=0; i
<natoms_solvt
; i
++,j
++)
361 combine_atoms(atoms
,atoms_solvt
,*x
,v
?*v
:NULL
,x_solvt
,v_solvt
,
362 &atoms_all
,&x_all
,&v_all
);
364 /* Do neighboursearching step */
365 do_nsgrid(stdout
,bVerbose
,box
,x_all
,atoms_all
,max_vdw
,oenv
);
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
++) {
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
);
380 for(j
=j0
; (j
<j1
&& nremove
<natoms_solvt
); j
++) {
381 jnr
= nlist
->jjnr
[j
];
382 copy_rvec(x_all
[jnr
],xj
);
384 /* Check solvent-protein and solvent-solvent */
385 is1
= inr
-natoms_prot
;
386 is2
= jnr
-natoms_prot
;
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
393 bXor((is1
>= 0),(is2
>= 0)) && /* One atom is protein */
394 ((is1
< 0) || ((is1
>= 0) && !remove
[is1
])) &&
395 ((is2
< 0) || ((is2
>= 0) && !remove
[is2
]))) ||
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 */
404 (atoms_solvt
->atom
[is1
].resind
!= /* Not the same residue */
405 atoms_solvt
->atom
[is2
].resind
))) {
410 r2
= sqr(r_all
[inr
]+r_all
[jnr
]);
413 nremove
= natoms_solvt
;
414 for(k
=0; k
<nremove
; k
++) {
418 /* Need only remove one of the solvents... */
420 (void) mark_res(is2
,remove
,natoms_solvt
,atoms_solvt
->atom
,
423 (void) mark_res(is1
,remove
,natoms_solvt
,atoms_solvt
->atom
,
426 fprintf(stderr
,"Neither atom is solvent%d %d\n",is1
,is2
);
432 fprintf(stderr
," tested %d pairs, removed %d atoms.\n",ntest
,nremove
);
436 for(i
=0; i
<natoms_solvt
; i
++)
437 fprintf(debug
,"remove[%5d] = %s\n",i
,bool_names
[remove
[i
]]);
439 /* Search again, now with another cut-off */
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
);
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];
451 for(j
=j0
; j
<j1
; j
++) {
452 jnr
= nlist
->jjnr
[j
];
454 /* Check solvent-protein and solvent-solvent */
455 is1
= inr
-natoms_prot
;
456 is2
= jnr
-natoms_prot
;
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
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
);
468 fprintf(stderr
,"Keeping %d solvent atoms after proximity check\n",
470 for (i
=0; i
<natoms_solvt
; i
++)
471 remove
[i
] = remove
[i
] || !keep
[i
];
474 /* count how many atoms and residues will be added and make space */
477 jnres
= atoms_solvt
->nres
;
481 for (i
=0; ((i
<atoms_solvt
->nr
) &&
482 ((max_sol
== 0) || (jnres
< max_sol
))); i
++) {
486 (atoms_solvt
->atom
[i
].resind
!= atoms_solvt
->atom
[i
-1].resind
))
492 fprintf(debug
,"Will add %d atoms in %d residues\n",j
,jnres
);
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
)
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
);
513 /* add the selected atoms_solvt to atoms */
515 resnr
= atoms
->resinfo
[atoms
->atom
[atoms
->nr
-1].resind
].nr
;
521 for (i
=0; i
<atoms_solvt
->nr
; i
++) {
524 atoms_solvt
->atom
[i
].resind
!= atoms_solvt
->atom
[prev
].resind
) {
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
]);
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;
547 srenew(atoms
->resinfo
, atoms
->nres
+nresadd
);
550 fprintf(stderr
,"Added %d molecules\n",nresadd
);
553 done_atom(atoms_all
);