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
,
287 TRUE
,FALSE
,FALSE
,NULL
);
290 dump_nblist(debug
,cr
,fr
,0);
293 fprintf(stderr
,"Successfully made neighbourlist\n");
296 gmx_bool
bXor(gmx_bool b1
,gmx_bool b2
)
298 return (b1
&& !b2
) || (b2
&& !b1
);
301 void add_conf(t_atoms
*atoms
, rvec
**x
, rvec
**v
, real
**r
, gmx_bool bSrenew
,
302 int ePBC
, matrix box
, gmx_bool bInsert
,
303 t_atoms
*atoms_solvt
,rvec
*x_solvt
,rvec
*v_solvt
,real
*r_solvt
,
304 gmx_bool bVerbose
,real rshell
,int max_sol
, const output_env_t oenv
)
308 real max_vdw
,*r_prot
,*r_all
,n2
,r2
,ib1
,ib2
;
309 int natoms_prot
,natoms_solvt
;
310 int i
,j
,jj
,m
,j0
,j1
,jjj
,jnres
,jnr
,inr
,iprot
,is1
,is2
;
311 int prev
,resnr
,nresadd
,d
,k
,ncells
,maxincell
;
312 int dx0
,dx1
,dy0
,dy1
,dz0
,dz1
;
313 int ntest
,nremove
,nkeep
;
314 rvec dx
,xi
,xj
,xpp
,*x_all
,*v_all
;
315 gmx_bool
*remove
,*keep
;
318 natoms_prot
= atoms
->nr
;
319 natoms_solvt
= atoms_solvt
->nr
;
320 if (natoms_solvt
<= 0) {
321 fprintf(stderr
,"WARNING: Nothing to add\n");
325 if (ePBC
== epbcSCREW
)
326 gmx_fatal(FARGS
,"Sorry, %s pbc is not yet supported",epbc_names
[ePBC
]);
329 fprintf(stderr
,"Calculating Overlap...\n");
331 /* Set margin around box edges to largest solvent dimension.
332 * The maximum distance between atoms in a solvent molecule should
333 * be calculated. At the moment a fudge factor of 3 is used.
336 box_margin
= 3*find_max_real(natoms_solvt
,r_solvt
);
337 max_vdw
= max(3*find_max_real(natoms_prot
,r_prot
),box_margin
);
338 fprintf(stderr
,"box_margin = %g\n",box_margin
);
340 snew(remove
,natoms_solvt
);
344 for(i
=0; i
<atoms_solvt
->nr
; i
++)
345 if ( outside_box_plus_margin(x_solvt
[i
],box
) )
346 i
=mark_res(i
,remove
,atoms_solvt
->nr
,atoms_solvt
->atom
,&nremove
);
347 fprintf(stderr
,"Removed %d atoms that were outside the box\n",nremove
);
350 /* Define grid stuff for genbox */
351 /* Largest VDW radius */
352 snew(r_all
,natoms_prot
+natoms_solvt
);
353 for(i
=j
=0; i
<natoms_prot
; i
++,j
++)
355 for(i
=0; i
<natoms_solvt
; i
++,j
++)
359 combine_atoms(atoms
,atoms_solvt
,*x
,v
?*v
:NULL
,x_solvt
,v_solvt
,
360 &atoms_all
,&x_all
,&v_all
);
362 /* Do neighboursearching step */
363 do_nsgrid(stdout
,bVerbose
,box
,x_all
,atoms_all
,max_vdw
,oenv
);
365 /* check solvent with solute */
366 nlist
= &(fr
->nblists
[0].nlist_sr
[eNL_VDW
]);
367 fprintf(stderr
,"nri = %d, nrj = %d\n",nlist
->nri
,nlist
->nrj
);
368 for(bSolSol
=0; (bSolSol
<=(bInsert
? 0 : 1)); bSolSol
++) {
370 fprintf(stderr
,"Checking %s-Solvent overlap:",
371 bSolSol
? "Solvent" : "Protein");
372 for(i
=0; (i
<nlist
->nri
&& nremove
<natoms_solvt
); i
++) {
373 inr
= nlist
->iinr
[i
];
374 j0
= nlist
->jindex
[i
];
375 j1
= nlist
->jindex
[i
+1];
376 rvec_add(x_all
[inr
],fr
->shift_vec
[nlist
->shift
[i
]],xi
);
378 for(j
=j0
; (j
<j1
&& nremove
<natoms_solvt
); j
++) {
379 jnr
= nlist
->jjnr
[j
];
380 copy_rvec(x_all
[jnr
],xj
);
382 /* Check solvent-protein and solvent-solvent */
383 is1
= inr
-natoms_prot
;
384 is2
= jnr
-natoms_prot
;
386 /* Check if at least one of the atoms is a solvent that is not yet
387 * listed for removal, and if both are solvent, that they are not in the
391 bXor((is1
>= 0),(is2
>= 0)) && /* One atom is protein */
392 ((is1
< 0) || ((is1
>= 0) && !remove
[is1
])) &&
393 ((is2
< 0) || ((is2
>= 0) && !remove
[is2
]))) ||
396 (is1
>= 0) && (!remove
[is1
]) && /* is1 is solvent */
397 (is2
>= 0) && (!remove
[is2
]) && /* is2 is solvent */
398 (bInsert
|| /* when inserting also check inside the box */
399 (outside_box_minus_margin2(x_solvt
[is1
],box
) && /* is1 on edge */
400 outside_box_minus_margin2(x_solvt
[is2
],box
)) /* is2 on edge */
402 (atoms_solvt
->atom
[is1
].resind
!= /* Not the same residue */
403 atoms_solvt
->atom
[is2
].resind
))) {
408 r2
= sqr(r_all
[inr
]+r_all
[jnr
]);
411 nremove
= natoms_solvt
;
412 for(k
=0; k
<nremove
; k
++) {
416 /* Need only remove one of the solvents... */
418 (void) mark_res(is2
,remove
,natoms_solvt
,atoms_solvt
->atom
,
421 (void) mark_res(is1
,remove
,natoms_solvt
,atoms_solvt
->atom
,
424 fprintf(stderr
,"Neither atom is solvent%d %d\n",is1
,is2
);
430 fprintf(stderr
," tested %d pairs, removed %d atoms.\n",ntest
,nremove
);
434 for(i
=0; i
<natoms_solvt
; i
++)
435 fprintf(debug
,"remove[%5d] = %s\n",i
,bool_names
[remove
[i
]]);
437 /* Search again, now with another cut-off */
439 do_nsgrid(stdout
,bVerbose
,box
,x_all
,atoms_all
,rshell
,oenv
);
440 nlist
= &(fr
->nblists
[0].nlist_sr
[eNL_VDW
]);
441 fprintf(stderr
,"nri = %d, nrj = %d\n",nlist
->nri
,nlist
->nrj
);
443 snew(keep
,natoms_solvt
);
444 for(i
=0; i
<nlist
->nri
; i
++) {
445 inr
= nlist
->iinr
[i
];
446 j0
= nlist
->jindex
[i
];
447 j1
= nlist
->jindex
[i
+1];
449 for(j
=j0
; j
<j1
; j
++) {
450 jnr
= nlist
->jjnr
[j
];
452 /* Check solvent-protein and solvent-solvent */
453 is1
= inr
-natoms_prot
;
454 is2
= jnr
-natoms_prot
;
456 /* Check if at least one of the atoms is a solvent that is not yet
457 * listed for removal, and if both are solvent, that they are not in the
461 mark_res(is1
,keep
,natoms_solvt
,atoms_solvt
->atom
,&nkeep
);
462 else if (is1
<0 && is2
>=0)
463 mark_res(is2
,keep
,natoms_solvt
,atoms_solvt
->atom
,&nkeep
);
466 fprintf(stderr
,"Keeping %d solvent atoms after proximity check\n",
468 for (i
=0; i
<natoms_solvt
; i
++)
469 remove
[i
] = remove
[i
] || !keep
[i
];
472 /* count how many atoms and residues will be added and make space */
475 jnres
= atoms_solvt
->nres
;
479 for (i
=0; ((i
<atoms_solvt
->nr
) &&
480 ((max_sol
== 0) || (jnres
< max_sol
))); i
++) {
484 (atoms_solvt
->atom
[i
].resind
!= atoms_solvt
->atom
[i
-1].resind
))
490 fprintf(debug
,"Will add %d atoms in %d residues\n",j
,jnres
);
492 /* Flag the remaing solvent atoms to be removed */
493 jjj
= atoms_solvt
->atom
[i
-1].resind
;
494 for ( ; (i
<atoms_solvt
->nr
); i
++) {
495 if (atoms_solvt
->atom
[i
].resind
> jjj
)
503 srenew(atoms
->resinfo
, atoms
->nres
+jnres
);
504 srenew(atoms
->atomname
, atoms
->nr
+j
);
505 srenew(atoms
->atom
, atoms
->nr
+j
);
506 srenew(*x
, atoms
->nr
+j
);
507 if (v
) srenew(*v
, atoms
->nr
+j
);
508 srenew(*r
, atoms
->nr
+j
);
511 /* add the selected atoms_solvt to atoms */
513 resnr
= atoms
->resinfo
[atoms
->atom
[atoms
->nr
-1].resind
].nr
;
519 for (i
=0; i
<atoms_solvt
->nr
; i
++) {
522 atoms_solvt
->atom
[i
].resind
!= atoms_solvt
->atom
[prev
].resind
) {
526 atoms
->resinfo
[atoms
->nres
-1] =
527 atoms_solvt
->resinfo
[atoms_solvt
->atom
[i
].resind
];
528 atoms
->resinfo
[atoms
->nres
-1].nr
= resnr
;
529 /* calculate shift of the solvent molecule using the first atom */
530 copy_rvec(x_solvt
[i
],dx
);
531 put_atoms_in_box(ePBC
,box
,1,&dx
);
532 rvec_dec(dx
,x_solvt
[i
]);
534 atoms
->atom
[atoms
->nr
] = atoms_solvt
->atom
[i
];
535 atoms
->atomname
[atoms
->nr
] = atoms_solvt
->atomname
[i
];
536 rvec_add(x_solvt
[i
],dx
,(*x
)[atoms
->nr
]);
537 if (v
) copy_rvec(v_solvt
[i
],(*v
)[atoms
->nr
]);
538 (*r
)[atoms
->nr
] = r_solvt
[i
];
539 atoms
->atom
[atoms
->nr
].resind
= atoms
->nres
-1;
545 srenew(atoms
->resinfo
, atoms
->nres
+nresadd
);
548 fprintf(stderr
,"Added %d molecules\n",nresadd
);
551 done_atom(atoms_all
);