?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Parallel vs. sequentiual code: I get binary identical result (apart from comment...
[gromacs/qmmm-gamess-us.git] / src / tools / addconf.c
blob6380cacc59799ae598334274beb50fe0f241c57b
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 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 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, 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,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,"Succesfully made neighbourlist\n");
311 }
312
313 bool bXor(bool b1,bool b2)
314 {
315 return (b1 && !b2) || (b2 && !b1);
316 }
317
318 void add_conf(t_atoms *atoms, rvec **x, rvec **v, real **r, bool bSrenew,
319 int ePBC, matrix box, bool bInsert,
320 t_atoms *atoms_solvt,rvec *x_solvt,rvec *v_solvt,real *r_solvt,
321 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 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 }
Atempt to write QMMM interface for Gamess-US