?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added Verlet scheme and NxN non-bonded functionality
[gromacs.git] / src / mdlib / constr.c
blobf59b3b69960753587587218d36d0a033653dbb1d
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROwing Monsters And Cloning Shrimps
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include "confio.h"
41 #include "constr.h"
42 #include "copyrite.h"
43 #include "invblock.h"
44 #include "main.h"
45 #include "mdrun.h"
46 #include "nrnb.h"
47 #include "smalloc.h"
48 #include "vec.h"
49 #include "physics.h"
50 #include "names.h"
51 #include "txtdump.h"
52 #include "domdec.h"
53 #include "pdbio.h"
54 #include "partdec.h"
55 #include "splitter.h"
56 #include "mtop_util.h"
57 #include "gmxfio.h"
58 #include "gmx_omp_nthreads.h"
59
60 typedef struct gmx_constr {
61 int ncon_tot; /* The total number of constraints */
62 int nflexcon; /* The number of flexible constraints */
63 int n_at2con_mt; /* The size of at2con = #moltypes */
64 t_blocka *at2con_mt; /* A list of atoms to constraints */
65 int n_at2settle_mt; /* The size of at2settle = #moltypes */
66 int **at2settle_mt; /* A list of atoms to settles */
67 gmx_bool bInterCGsettles;
68 gmx_lincsdata_t lincsd; /* LINCS data */
69 gmx_shakedata_t shaked; /* SHAKE data */
70 gmx_settledata_t settled; /* SETTLE data */
71 int nblocks; /* The number of SHAKE blocks */
72 int *sblock; /* The SHAKE blocks */
73 int sblock_nalloc;/* The allocation size of sblock */
74 real *lagr; /* Lagrange multipliers for SHAKE */
75 int lagr_nalloc; /* The allocation size of lagr */
76 int maxwarn; /* The maximum number of warnings */
77 int warncount_lincs;
78 int warncount_settle;
79 gmx_edsam_t ed; /* The essential dynamics data */
80
81 tensor *rmdr_th; /* Thread local working data */
82 int *settle_error; /* Thread local working data */
83
84 gmx_mtop_t *warn_mtop; /* Only used for printing warnings */
85 } t_gmx_constr;
86
87 typedef struct {
88 atom_id iatom[3];
89 atom_id blocknr;
90 } t_sortblock;
91
92 static void *init_vetavars(t_vetavars *vars,
93 gmx_bool constr_deriv,
94 real veta,real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
95 {
96 double g;
97 int i;
98
99 /* first, set the alpha integrator variable */
100 if ((ir->opts.nrdf[0] > 0) && bPscal)
101 {
102 vars->alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
103 } else {
104 vars->alpha = 1.0;
105 }
106 g = 0.5*veta*ir->delta_t;
107 vars->rscale = exp(g)*series_sinhx(g);
108 g = -0.25*vars->alpha*veta*ir->delta_t;
109 vars->vscale = exp(g)*series_sinhx(g);
110 vars->rvscale = vars->vscale*vars->rscale;
111 vars->veta = vetanew;
112
113 if (constr_deriv)
114 {
115 snew(vars->vscale_nhc,ir->opts.ngtc);
116 if ((ekind==NULL) || (!bPscal))
117 {
118 for (i=0;i<ir->opts.ngtc;i++)
119 {
120 vars->vscale_nhc[i] = 1;
121 }
122 }
123 else
124 {
125 for (i=0;i<ir->opts.ngtc;i++)
126 {
127 vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
128 }
129 }
130 }
131 else
132 {
133 vars->vscale_nhc = NULL;
134 }
135
136 return vars;
137 }
138
139 static void free_vetavars(t_vetavars *vars)
140 {
141 if (vars->vscale_nhc != NULL)
142 {
143 sfree(vars->vscale_nhc);
144 }
145 }
146
147 static int pcomp(const void *p1, const void *p2)
148 {
149 int db;
150 atom_id min1,min2,max1,max2;
151 t_sortblock *a1=(t_sortblock *)p1;
152 t_sortblock *a2=(t_sortblock *)p2;
153
154 db=a1->blocknr-a2->blocknr;
155
156 if (db != 0)
157 return db;
158
159 min1=min(a1->iatom[1],a1->iatom[2]);
160 max1=max(a1->iatom[1],a1->iatom[2]);
161 min2=min(a2->iatom[1],a2->iatom[2]);
162 max2=max(a2->iatom[1],a2->iatom[2]);
163
164 if (min1 == min2)
165 return max1-max2;
166 else
167 return min1-min2;
168 }
169
170 int n_flexible_constraints(struct gmx_constr *constr)
171 {
172 int nflexcon;
173
174 if (constr)
175 nflexcon = constr->nflexcon;
176 else
177 nflexcon = 0;
178
179 return nflexcon;
180 }
181
182 void too_many_constraint_warnings(int eConstrAlg,int warncount)
183 {
184 const char *abort="- aborting to avoid logfile runaway.\n"
185 "This normally happens when your system is not sufficiently equilibrated,"
186 "or if you are changing lambda too fast in free energy simulations.\n";
187
188 gmx_fatal(FARGS,
189 "Too many %s warnings (%d)\n"
190 "If you know what you are doing you can %s"
191 "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
192 "but normally it is better to fix the problem",
193 (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE",warncount,
194 (eConstrAlg == econtLINCS) ?
195 "adjust the lincs warning threshold in your mdp file\nor " : "\n");
196 }
197
198 static void write_constr_pdb(const char *fn,const char *title,
199 gmx_mtop_t *mtop,
200 int start,int homenr,t_commrec *cr,
201 rvec x[],matrix box)
202 {
203 char fname[STRLEN],format[STRLEN];
204 FILE *out;
205 int dd_ac0=0,dd_ac1=0,i,ii,resnr;
206 gmx_domdec_t *dd;
207 char *anm,*resnm;
208
209 dd = NULL;
210 if (PAR(cr))
211 {
212 sprintf(fname,"%s_n%d.pdb",fn,cr->sim_nodeid);
213 if (DOMAINDECOMP(cr))
214 {
215 dd = cr->dd;
216 dd_get_constraint_range(dd,&dd_ac0,&dd_ac1);
217 start = 0;
218 homenr = dd_ac1;
219 }
220 }
221 else
222 {
223 sprintf(fname,"%s.pdb",fn);
224 }
225 sprintf(format,"%s\n",pdbformat);
226
227 out = gmx_fio_fopen(fname,"w");
228
229 fprintf(out,"TITLE %s\n",title);
230 gmx_write_pdb_box(out,-1,box);
231 for(i=start; i<start+homenr; i++)
232 {
233 if (dd != NULL)
234 {
235 if (i >= dd->nat_home && i < dd_ac0)
236 {
237 continue;
238 }
239 ii = dd->gatindex[i];
240 }
241 else
242 {
243 ii = i;
244 }
245 gmx_mtop_atominfo_global(mtop,ii,&anm,&resnr,&resnm);
246 fprintf(out,format,"ATOM",(ii+1)%100000,
247 anm,resnm,' ',resnr%10000,' ',
248 10*x[i][XX],10*x[i][YY],10*x[i][ZZ]);
249 }
250 fprintf(out,"TER\n");
251
252 gmx_fio_fclose(out);
253 }
254
255 static void dump_confs(FILE *fplog,gmx_large_int_t step,gmx_mtop_t *mtop,
256 int start,int homenr,t_commrec *cr,
257 rvec x[],rvec xprime[],matrix box)
258 {
259 char buf[256],buf2[22];
260
261 char *env=getenv("GMX_SUPPRESS_DUMP");
262 if (env)
263 return;
264
265 sprintf(buf,"step%sb",gmx_step_str(step,buf2));
266 write_constr_pdb(buf,"initial coordinates",
267 mtop,start,homenr,cr,x,box);
268 sprintf(buf,"step%sc",gmx_step_str(step,buf2));
269 write_constr_pdb(buf,"coordinates after constraining",
270 mtop,start,homenr,cr,xprime,box);
271 if (fplog)
272 {
273 fprintf(fplog,"Wrote pdb files with previous and current coordinates\n");
274 }
275 fprintf(stderr,"Wrote pdb files with previous and current coordinates\n");
276 }
277
278 static void pr_sortblock(FILE *fp,const char *title,int nsb,t_sortblock sb[])
279 {
280 int i;
281
282 fprintf(fp,"%s\n",title);
283 for(i=0; (i<nsb); i++)
284 fprintf(fp,"i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
285 i,sb[i].iatom[0],sb[i].iatom[1],sb[i].iatom[2],
286 sb[i].blocknr);
287 }
288
289 gmx_bool constrain(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
290 struct gmx_constr *constr,
291 t_idef *idef,t_inputrec *ir,gmx_ekindata_t *ekind,
292 t_commrec *cr,
293 gmx_large_int_t step,int delta_step,
294 t_mdatoms *md,
295 rvec *x,rvec *xprime,rvec *min_proj,
296 gmx_bool bMolPBC,matrix box,
297 real lambda,real *dvdlambda,
298 rvec *v,tensor *vir,
299 t_nrnb *nrnb,int econq,gmx_bool bPscal,
300 real veta, real vetanew)
301 {
302 gmx_bool bOK,bDump;
303 int start,homenr,nrend;
304 int i,j,d;
305 int ncons,settle_error;
306 tensor rmdr;
307 rvec *vstor;
308 real invdt,vir_fac,t;
309 t_ilist *settle;
310 int nsettle;
311 t_pbc pbc,*pbc_null;
312 char buf[22];
313 t_vetavars vetavar;
314 int nth,th;
315
316 if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
317 {
318 gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
319 }
320
321 bOK = TRUE;
322 bDump = FALSE;
323
324 start = md->start;
325 homenr = md->homenr;
326 nrend = start+homenr;
327
328 /* set constants for pressure control integration */
329 init_vetavars(&vetavar,econq!=econqCoord,
330 veta,vetanew,ir,ekind,bPscal);
331
332 if (ir->delta_t == 0)
333 {
334 invdt = 0;
335 }
336 else
337 {
338 invdt = 1/ir->delta_t;
339 }
340
341 if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
342 {
343 /* Set the constraint lengths for the step at which this configuration
344 * is meant to be. The invmasses should not be changed.
345 */
346 lambda += delta_step*ir->fepvals->delta_lambda;
347 }
348
349 if (vir != NULL)
350 {
351 clear_mat(rmdr);
352 }
353
354 where();
355
356 settle = &idef->il[F_SETTLE];
357 nsettle = settle->nr/(1+NRAL(F_SETTLE));
358
359 if (nsettle > 0)
360 {
361 nth = gmx_omp_nthreads_get(emntSETTLE);
362 }
363 else
364 {
365 nth = 1;
366 }
367
368 if (nth > 1 && constr->rmdr_th == NULL)
369 {
370 snew(constr->rmdr_th,nth);
371 snew(constr->settle_error,nth);
372 }
373
374 settle_error = -1;
375
376 /* We do not need full pbc when constraints do not cross charge groups,
377 * i.e. when dd->constraint_comm==NULL.
378 * Note that PBC for constraints is different from PBC for bondeds.
379 * For constraints there is both forward and backward communication.
380 */
381 if (ir->ePBC != epbcNONE &&
382 (cr->dd || bMolPBC) && !(cr->dd && cr->dd->constraint_comm==NULL))
383 {
384 /* With pbc=screw the screw has been changed to a shift
385 * by the constraint coordinate communication routine,
386 * so that here we can use normal pbc.
387 */
388 pbc_null = set_pbc_dd(&pbc,ir->ePBC,cr->dd,FALSE,box);
389 }
390 else
391 {
392 pbc_null = NULL;
393 }
394
395 /* Communicate the coordinates required for the non-local constraints
396 * for LINCS and/or SETTLE.
397 */
398 if (cr->dd)
399 {
400 dd_move_x_constraints(cr->dd,box,x,xprime);
401 }
402 else if (PARTDECOMP(cr))
403 {
404 pd_move_x_constraints(cr,x,xprime);
405 }
406
407 if (constr->lincsd != NULL)
408 {
409 bOK = constrain_lincs(fplog,bLog,bEner,ir,step,constr->lincsd,md,cr,
410 x,xprime,min_proj,
411 box,pbc_null,lambda,dvdlambda,
412 invdt,v,vir!=NULL,rmdr,
413 econq,nrnb,
414 constr->maxwarn,&constr->warncount_lincs);
415 if (!bOK && constr->maxwarn >= 0)
416 {
417 if (fplog != NULL)
418 {
419 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
420 econstr_names[econtLINCS],gmx_step_str(step,buf));
421 }
422 bDump = TRUE;
423 }
424 }
425
426 if (constr->nblocks > 0)
427 {
428 switch (econq) {
429 case (econqCoord):
430 bOK = bshakef(fplog,constr->shaked,
431 homenr,md->invmass,constr->nblocks,constr->sblock,
432 idef,ir,x,xprime,nrnb,
433 constr->lagr,lambda,dvdlambda,
434 invdt,v,vir!=NULL,rmdr,constr->maxwarn>=0,econq,&vetavar);
435 break;
436 case (econqVeloc):
437 bOK = bshakef(fplog,constr->shaked,
438 homenr,md->invmass,constr->nblocks,constr->sblock,
439 idef,ir,x,min_proj,nrnb,
440 constr->lagr,lambda,dvdlambda,
441 invdt,NULL,vir!=NULL,rmdr,constr->maxwarn>=0,econq,&vetavar);
442 break;
443 default:
444 gmx_fatal(FARGS,"Internal error, SHAKE called for constraining something else than coordinates");
445 break;
446 }
447
448 if (!bOK && constr->maxwarn >= 0)
449 {
450 if (fplog != NULL)
451 {
452 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
453 econstr_names[econtSHAKE],gmx_step_str(step,buf));
454 }
455 bDump = TRUE;
456 }
457 }
458
459 if (nsettle > 0)
460 {
461 int calcvir_atom_end;
462
463 if (vir == NULL)
464 {
465 calcvir_atom_end = 0;
466 }
467 else
468 {
469 calcvir_atom_end = md->start + md->homenr;
470 }
471
472 switch (econq)
473 {
474 case econqCoord:
475 #pragma omp parallel for num_threads(nth) schedule(static)
476 for(th=0; th<nth; th++)
477 {
478 int start_th,end_th;
479
480 if (th > 0)
481 {
482 clear_mat(constr->rmdr_th[th]);
483 }
484
485 start_th = (nsettle* th )/nth;
486 end_th = (nsettle*(th+1))/nth;
487 if (start_th >= 0 && end_th - start_th > 0)
488 {
489 csettle(constr->settled,
490 end_th-start_th,
491 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
492 pbc_null,
493 x[0],xprime[0],
494 invdt,v?v[0]:NULL,calcvir_atom_end,
495 th == 0 ? rmdr : constr->rmdr_th[th],
496 th == 0 ? &settle_error : &constr->settle_error[th],
497 &vetavar);
498 }
499 }
500 inc_nrnb(nrnb,eNR_SETTLE,nsettle);
501 if (v != NULL)
502 {
503 inc_nrnb(nrnb,eNR_CONSTR_V,nsettle*3);
504 }
505 if (vir != NULL)
506 {
507 inc_nrnb(nrnb,eNR_CONSTR_VIR,nsettle*3);
508 }
509 break;
510 case econqVeloc:
511 case econqDeriv:
512 case econqForce:
513 case econqForceDispl:
514 #pragma omp parallel for num_threads(nth) schedule(static)
515 for(th=0; th<nth; th++)
516 {
517 int start_th,end_th;
518
519 if (th > 0)
520 {
521 clear_mat(constr->rmdr_th[th]);
522 }
523
524 start_th = (nsettle* th )/nth;
525 end_th = (nsettle*(th+1))/nth;
526
527 if (start_th >= 0 && end_th - start_th > 0)
528 {
529 settle_proj(fplog,constr->settled,econq,
530 end_th-start_th,
531 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
532 pbc_null,
533 x,
534 xprime,min_proj,calcvir_atom_end,
535 th == 0 ? rmdr : constr->rmdr_th[th],
536 &vetavar);
537 }
538 }
539 /* This is an overestimate */
540 inc_nrnb(nrnb,eNR_SETTLE,nsettle);
541 break;
542 case econqDeriv_FlexCon:
543 /* Nothing to do, since the are no flexible constraints in settles */
544 break;
545 default:
546 gmx_incons("Unknown constraint quantity for settle");
547 }
548 }
549
550 if (settle->nr > 0)
551 {
552 /* Combine virial and error info of the other threads */
553 for(i=1; i<nth; i++)
554 {
555 m_add(rmdr,constr->rmdr_th[i],rmdr);
556 settle_error = constr->settle_error[i];
557 }
558
559 if (econq == econqCoord && settle_error >= 0)
560 {
561 bOK = FALSE;
562 if (constr->maxwarn >= 0)
563 {
564 char buf[256];
565 sprintf(buf,
566 "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
567 "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
568 step,ddglatnr(cr->dd,settle->iatoms[settle_error*(1+NRAL(F_SETTLE))+1]));
569 if (fplog)
570 {
571 fprintf(fplog,"%s",buf);
572 }
573 fprintf(stderr,"%s",buf);
574 constr->warncount_settle++;
575 if (constr->warncount_settle > constr->maxwarn)
576 {
577 too_many_constraint_warnings(-1,constr->warncount_settle);
578 }
579 bDump = TRUE;
580 }
581 }
582 }
583
584 free_vetavars(&vetavar);
585
586 if (vir != NULL)
587 {
588 switch (econq)
589 {
590 case econqCoord:
591 vir_fac = 0.5/(ir->delta_t*ir->delta_t);
592 break;
593 case econqVeloc:
594 vir_fac = 0.5/ir->delta_t;
595 break;
596 case econqForce:
597 case econqForceDispl:
598 vir_fac = 0.5;
599 break;
600 default:
601 vir_fac = 0;
602 gmx_incons("Unsupported constraint quantity for virial");
603 }
604
605 if (EI_VV(ir->eI))
606 {
607 vir_fac *= 2; /* only constraining over half the distance here */
608 }
609 for(i=0; i<DIM; i++)
610 {
611 for(j=0; j<DIM; j++)
612 {
613 (*vir)[i][j] = vir_fac*rmdr[i][j];
614 }
615 }
616 }
617
618 if (bDump)
619 {
620 dump_confs(fplog,step,constr->warn_mtop,start,homenr,cr,x,xprime,box);
621 }
622
623 if (econq == econqCoord)
624 {
625 if (ir->ePull == epullCONSTRAINT)
626 {
627 if (EI_DYNAMICS(ir->eI))
628 {
629 t = ir->init_t + (step + delta_step)*ir->delta_t;
630 }
631 else
632 {
633 t = ir->init_t;
634 }
635 set_pbc(&pbc,ir->ePBC,box);
636 pull_constraint(ir->pull,md,&pbc,cr,ir->delta_t,t,x,xprime,v,*vir);
637 }
638 if (constr->ed && delta_step > 0)
639 {
640 /* apply the essential dynamcs constraints here */
641 do_edsam(ir,step,md,cr,xprime,v,box,constr->ed);
642 }
643 }
644
645 return bOK;
646 }
647
648 real *constr_rmsd_data(struct gmx_constr *constr)
649 {
650 if (constr->lincsd)
651 return lincs_rmsd_data(constr->lincsd);
652 else
653 return NULL;
654 }
655
656 real constr_rmsd(struct gmx_constr *constr,gmx_bool bSD2)
657 {
658 if (constr->lincsd)
659 return lincs_rmsd(constr->lincsd,bSD2);
660 else
661 return 0;
662 }
663
664 static void make_shake_sblock_pd(struct gmx_constr *constr,
665 t_idef *idef,t_mdatoms *md)
666 {
667 int i,j,m,ncons;
668 int bstart,bnr;
669 t_blocka sblocks;
670 t_sortblock *sb;
671 t_iatom *iatom;
672 atom_id *inv_sblock;
673
674 /* Since we are processing the local topology,
675 * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
676 */
677 ncons = idef->il[F_CONSTR].nr/3;
678
679 init_blocka(&sblocks);
680 gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
681
682 /*
683 bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
684 nblocks=blocks->multinr[idef->nodeid] - bstart;
685 */
686 bstart = 0;
687 constr->nblocks = sblocks.nr;
688 if (debug)
689 fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
690 ncons,bstart,constr->nblocks);
691
692 /* Calculate block number for each atom */
693 inv_sblock = make_invblocka(&sblocks,md->nr);
694
695 done_blocka(&sblocks);
696
697 /* Store the block number in temp array and
698 * sort the constraints in order of the sblock number
699 * and the atom numbers, really sorting a segment of the array!
700 */
701 #ifdef DEBUGIDEF
702 pr_idef(fplog,0,"Before Sort",idef);
703 #endif
704 iatom=idef->il[F_CONSTR].iatoms;
705 snew(sb,ncons);
706 for(i=0; (i<ncons); i++,iatom+=3) {
707 for(m=0; (m<3); m++)
708 sb[i].iatom[m] = iatom[m];
709 sb[i].blocknr = inv_sblock[iatom[1]];
710 }
711
712 /* Now sort the blocks */
713 if (debug) {
714 pr_sortblock(debug,"Before sorting",ncons,sb);
715 fprintf(debug,"Going to sort constraints\n");
716 }
717
718 qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
719
720 if (debug) {
721 pr_sortblock(debug,"After sorting",ncons,sb);
722 }
723
724 iatom=idef->il[F_CONSTR].iatoms;
725 for(i=0; (i<ncons); i++,iatom+=3)
726 for(m=0; (m<3); m++)
727 iatom[m]=sb[i].iatom[m];
728 #ifdef DEBUGIDEF
729 pr_idef(fplog,0,"After Sort",idef);
730 #endif
731
732 j=0;
733 snew(constr->sblock,constr->nblocks+1);
734 bnr=-2;
735 for(i=0; (i<ncons); i++) {
736 if (sb[i].blocknr != bnr) {
737 bnr=sb[i].blocknr;
738 constr->sblock[j++]=3*i;
739 }
740 }
741 /* Last block... */
742 constr->sblock[j++] = 3*ncons;
743
744 if (j != (constr->nblocks+1)) {
745 fprintf(stderr,"bstart: %d\n",bstart);
746 fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
747 j,constr->nblocks,ncons);
748 for(i=0; (i<ncons); i++)
749 fprintf(stderr,"i: %5d sb[i].blocknr: %5u\n",i,sb[i].blocknr);
750 for(j=0; (j<=constr->nblocks); j++)
751 fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
752 gmx_fatal(FARGS,"DEATH HORROR: "
753 "sblocks does not match idef->il[F_CONSTR]");
754 }
755 sfree(sb);
756 sfree(inv_sblock);
757 }
758
759 static void make_shake_sblock_dd(struct gmx_constr *constr,
760 t_ilist *ilcon,t_block *cgs,
761 gmx_domdec_t *dd)
762 {
763 int ncons,c,cg;
764 t_iatom *iatom;
765
766 if (dd->ncg_home+1 > constr->sblock_nalloc) {
767 constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
768 srenew(constr->sblock,constr->sblock_nalloc);
769 }
770
771 ncons = ilcon->nr/3;
772 iatom = ilcon->iatoms;
773 constr->nblocks = 0;
774 cg = 0;
775 for(c=0; c<ncons; c++) {
776 if (c == 0 || iatom[1] >= cgs->index[cg+1]) {
777 constr->sblock[constr->nblocks++] = 3*c;
778 while (iatom[1] >= cgs->index[cg+1])
779 cg++;
780 }
781 iatom += 3;
782 }
783 constr->sblock[constr->nblocks] = 3*ncons;
784 }
785
786 t_blocka make_at2con(int start,int natoms,
787 t_ilist *ilist,t_iparams *iparams,
788 gmx_bool bDynamics,int *nflexiblecons)
789 {
790 int *count,ncon,con,con_tot,nflexcon,ftype,i,a;
791 t_iatom *ia;
792 t_blocka at2con;
793 gmx_bool bFlexCon;
794
795 snew(count,natoms);
796 nflexcon = 0;
797 for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
798 ncon = ilist[ftype].nr/3;
799 ia = ilist[ftype].iatoms;
800 for(con=0; con<ncon; con++) {
801 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
802 iparams[ia[0]].constr.dB == 0);
803 if (bFlexCon)
804 nflexcon++;
805 if (bDynamics || !bFlexCon) {
806 for(i=1; i<3; i++) {
807 a = ia[i] - start;
808 count[a]++;
809 }
810 }
811 ia += 3;
812 }
813 }
814 *nflexiblecons = nflexcon;
815
816 at2con.nr = natoms;
817 at2con.nalloc_index = at2con.nr+1;
818 snew(at2con.index,at2con.nalloc_index);
819 at2con.index[0] = 0;
820 for(a=0; a<natoms; a++) {
821 at2con.index[a+1] = at2con.index[a] + count[a];
822 count[a] = 0;
823 }
824 at2con.nra = at2con.index[natoms];
825 at2con.nalloc_a = at2con.nra;
826 snew(at2con.a,at2con.nalloc_a);
827
828 /* The F_CONSTRNC constraints have constraint numbers
829 * that continue after the last F_CONSTR constraint.
830 */
831 con_tot = 0;
832 for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
833 ncon = ilist[ftype].nr/3;
834 ia = ilist[ftype].iatoms;
835 for(con=0; con<ncon; con++) {
836 bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
837 iparams[ia[0]].constr.dB == 0);
838 if (bDynamics || !bFlexCon) {
839 for(i=1; i<3; i++) {
840 a = ia[i] - start;
841 at2con.a[at2con.index[a]+count[a]++] = con_tot;
842 }
843 }
844 con_tot++;
845 ia += 3;
846 }
847 }
848
849 sfree(count);
850
851 return at2con;
852 }
853
854 static int *make_at2settle(int natoms,const t_ilist *ilist)
855 {
856 int *at2s;
857 int a,stride,s;
858
859 snew(at2s,natoms);
860 /* Set all to no settle */
861 for(a=0; a<natoms; a++)
862 {
863 at2s[a] = -1;
864 }
865
866 stride = 1 + NRAL(F_SETTLE);
867
868 for(s=0; s<ilist->nr; s+=stride)
869 {
870 at2s[ilist->iatoms[s+1]] = s/stride;
871 at2s[ilist->iatoms[s+2]] = s/stride;
872 at2s[ilist->iatoms[s+3]] = s/stride;
873 }
874
875 return at2s;
876 }
877
878 void set_constraints(struct gmx_constr *constr,
879 gmx_localtop_t *top,t_inputrec *ir,
880 t_mdatoms *md,t_commrec *cr)
881 {
882 t_idef *idef;
883 int ncons;
884 t_ilist *settle;
885 int iO,iH;
886
887 idef = &top->idef;
888
889 if (constr->ncon_tot > 0)
890 {
891 /* We are using the local topology,
892 * so there are only F_CONSTR constraints.
893 */
894 ncons = idef->il[F_CONSTR].nr/3;
895
896 /* With DD we might also need to call LINCS with ncons=0 for
897 * communicating coordinates to other nodes that do have constraints.
898 */
899 if (ir->eConstrAlg == econtLINCS)
900 {
901 set_lincs(idef,md,EI_DYNAMICS(ir->eI),cr,constr->lincsd);
902 }
903 if (ir->eConstrAlg == econtSHAKE)
904 {
905 if (cr->dd)
906 {
907 make_shake_sblock_dd(constr,&idef->il[F_CONSTR],&top->cgs,cr->dd);
908 }
909 else
910 {
911 make_shake_sblock_pd(constr,idef,md);
912 }
913 if (ncons > constr->lagr_nalloc)
914 {
915 constr->lagr_nalloc = over_alloc_dd(ncons);
916 srenew(constr->lagr,constr->lagr_nalloc);
917 }
918 }
919 }
920
921 if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
922 {
923 settle = &idef->il[F_SETTLE];
924 iO = settle->iatoms[1];
925 iH = settle->iatoms[2];
926 constr->settled =
927 settle_init(md->massT[iO],md->massT[iH],
928 md->invmass[iO],md->invmass[iH],
929 idef->iparams[settle->iatoms[0]].settle.doh,
930 idef->iparams[settle->iatoms[0]].settle.dhh);
931 }
932
933 /* Make a selection of the local atoms for essential dynamics */
934 if (constr->ed && cr->dd)
935 {
936 dd_make_local_ed_indices(cr->dd,constr->ed);
937 }
938 }
939
940 static void constr_recur(t_blocka *at2con,
941 t_ilist *ilist,t_iparams *iparams,gmx_bool bTopB,
942 int at,int depth,int nc,int *path,
943 real r0,real r1,real *r2max,
944 int *count)
945 {
946 int ncon1;
947 t_iatom *ia1,*ia2;
948 int c,con,a1;
949 gmx_bool bUse;
950 t_iatom *ia;
951 real len,rn0,rn1;
952
953 (*count)++;
954
955 ncon1 = ilist[F_CONSTR].nr/3;
956 ia1 = ilist[F_CONSTR].iatoms;
957 ia2 = ilist[F_CONSTRNC].iatoms;
958
959 /* Loop over all constraints connected to this atom */
960 for(c=at2con->index[at]; c<at2con->index[at+1]; c++) {
961 con = at2con->a[c];
962 /* Do not walk over already used constraints */
963 bUse = TRUE;
964 for(a1=0; a1<depth; a1++) {
965 if (con == path[a1])
966 bUse = FALSE;
967 }
968 if (bUse) {
969 ia = constr_iatomptr(ncon1,ia1,ia2,con);
970 /* Flexible constraints currently have length 0, which is incorrect */
971 if (!bTopB)
972 len = iparams[ia[0]].constr.dA;
973 else
974 len = iparams[ia[0]].constr.dB;
975 /* In the worst case the bond directions alternate */
976 if (nc % 2 == 0) {
977 rn0 = r0 + len;
978 rn1 = r1;
979 } else {
980 rn0 = r0;
981 rn1 = r1 + len;
982 }
983 /* Assume angles of 120 degrees between all bonds */
984 if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max) {
985 *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
986 if (debug) {
987 fprintf(debug,"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,rn1,sqrt(*r2max));
988 for(a1=0; a1<depth; a1++)
989 fprintf(debug," %d %5.3f",
990 path[a1],
991 iparams[constr_iatomptr(ncon1,ia1,ia2,con)[0]].constr.dA);
992 fprintf(debug," %d %5.3f\n",con,len);
993 }
994 }
995 /* Limit the number of recursions to 1000*nc,
996 * so a call does not take more than a second,
997 * even for highly connected systems.
998 */
999 if (depth + 1 < nc && *count < 1000*nc) {
1000 if (ia[1] == at)
1001 a1 = ia[2];
1002 else
1003 a1 = ia[1];
1004 /* Recursion */
1005 path[depth] = con;
1006 constr_recur(at2con,ilist,iparams,
1007 bTopB,a1,depth+1,nc,path,rn0,rn1,r2max,count);
1008 path[depth] = -1;
1009 }
1010 }
1011 }
1012 }
1013
1014 static real constr_r_max_moltype(FILE *fplog,
1015 gmx_moltype_t *molt,t_iparams *iparams,
1016 t_inputrec *ir)
1017 {
1018 int natoms,nflexcon,*path,at,count;
1019
1020 t_blocka at2con;
1021 real r0,r1,r2maxA,r2maxB,rmax,lam0,lam1;
1022
1023 if (molt->ilist[F_CONSTR].nr == 0 &&
1024 molt->ilist[F_CONSTRNC].nr == 0) {
1025 return 0;
1026 }
1027
1028 natoms = molt->atoms.nr;
1029
1030 at2con = make_at2con(0,natoms,molt->ilist,iparams,
1031 EI_DYNAMICS(ir->eI),&nflexcon);
1032 snew(path,1+ir->nProjOrder);
1033 for(at=0; at<1+ir->nProjOrder; at++)
1034 path[at] = -1;
1035
1036 r2maxA = 0;
1037 for(at=0; at<natoms; at++) {
1038 r0 = 0;
1039 r1 = 0;
1040
1041 count = 0;
1042 constr_recur(&at2con,molt->ilist,iparams,
1043 FALSE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxA,&count);
1044 }
1045 if (ir->efep == efepNO) {
1046 rmax = sqrt(r2maxA);
1047 } else {
1048 r2maxB = 0;
1049 for(at=0; at<natoms; at++) {
1050 r0 = 0;
1051 r1 = 0;
1052 count = 0;
1053 constr_recur(&at2con,molt->ilist,iparams,
1054 TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
1055 }
1056 lam0 = ir->fepvals->init_lambda;
1057 if (EI_DYNAMICS(ir->eI))
1058 lam0 += ir->init_step*ir->fepvals->delta_lambda;
1059 rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
1060 if (EI_DYNAMICS(ir->eI)) {
1061 lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
1062 rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
1063 }
1064 }
1065
1066 done_blocka(&at2con);
1067 sfree(path);
1068
1069 return rmax;
1070 }
1071
1072 real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir)
1073 {
1074 int mt;
1075 real rmax;
1076
1077 rmax = 0;
1078 for(mt=0; mt<mtop->nmoltype; mt++) {
1079 rmax = max(rmax,
1080 constr_r_max_moltype(fplog,&mtop->moltype[mt],
1081 mtop->ffparams.iparams,ir));
1082 }
1083
1084 if (fplog)
1085 fprintf(fplog,"Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n",1+ir->nProjOrder,rmax);
1086
1087 return rmax;
1088 }
1089
1090 gmx_constr_t init_constraints(FILE *fplog,
1091 gmx_mtop_t *mtop,t_inputrec *ir,
1092 gmx_edsam_t ed,t_state *state,
1093 t_commrec *cr)
1094 {
1095 int ncon,nset,nmol,settle_type,i,natoms,mt,nflexcon;
1096 struct gmx_constr *constr;
1097 char *env;
1098 t_ilist *ilist;
1099 gmx_mtop_ilistloop_t iloop;
1100
1101 ncon =
1102 gmx_mtop_ftype_count(mtop,F_CONSTR) +
1103 gmx_mtop_ftype_count(mtop,F_CONSTRNC);
1104 nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
1105
1106 if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
1107 {
1108 return NULL;
1109 }
1110
1111 snew(constr,1);
1112
1113 constr->ncon_tot = ncon;
1114 constr->nflexcon = 0;
1115 if (ncon > 0)
1116 {
1117 constr->n_at2con_mt = mtop->nmoltype;
1118 snew(constr->at2con_mt,constr->n_at2con_mt);
1119 for(mt=0; mt<mtop->nmoltype; mt++)
1120 {
1121 constr->at2con_mt[mt] = make_at2con(0,mtop->moltype[mt].atoms.nr,
1122 mtop->moltype[mt].ilist,
1123 mtop->ffparams.iparams,
1124 EI_DYNAMICS(ir->eI),&nflexcon);
1125 for(i=0; i<mtop->nmolblock; i++)
1126 {
1127 if (mtop->molblock[i].type == mt)
1128 {
1129 constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
1130 }
1131 }
1132 }
1133
1134 if (constr->nflexcon > 0)
1135 {
1136 if (fplog)
1137 {
1138 fprintf(fplog,"There are %d flexible constraints\n",
1139 constr->nflexcon);
1140 if (ir->fc_stepsize == 0)
1141 {
1142 fprintf(fplog,"\n"
1143 "WARNING: step size for flexible constraining = 0\n"
1144 " All flexible constraints will be rigid.\n"
1145 " Will try to keep all flexible constraints at their original length,\n"
1146 " but the lengths may exhibit some drift.\n\n");
1147 constr->nflexcon = 0;
1148 }
1149 }
1150 if (constr->nflexcon > 0)
1151 {
1152 please_cite(fplog,"Hess2002");
1153 }
1154 }
1155
1156 if (ir->eConstrAlg == econtLINCS)
1157 {
1158 constr->lincsd = init_lincs(fplog,mtop,
1159 constr->nflexcon,constr->at2con_mt,
1160 DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1161 ir->nLincsIter,ir->nProjOrder);
1162 }
1163
1164 if (ir->eConstrAlg == econtSHAKE) {
1165 if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1166 {
1167 gmx_fatal(FARGS,"SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1168 }
1169 if (constr->nflexcon)
1170 {
1171 gmx_fatal(FARGS,"For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1172 }
1173 please_cite(fplog,"Ryckaert77a");
1174 if (ir->bShakeSOR)
1175 {
1176 please_cite(fplog,"Barth95a");
1177 }
1178
1179 constr->shaked = shake_init();
1180 }
1181 }
1182
1183 if (nset > 0) {
1184 please_cite(fplog,"Miyamoto92a");
1185
1186 constr->bInterCGsettles = inter_charge_group_settles(mtop);
1187
1188 /* Check that we have only one settle type */
1189 settle_type = -1;
1190 iloop = gmx_mtop_ilistloop_init(mtop);
1191 while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
1192 {
1193 for (i=0; i<ilist[F_SETTLE].nr; i+=4)
1194 {
1195 if (settle_type == -1)
1196 {
1197 settle_type = ilist[F_SETTLE].iatoms[i];
1198 }
1199 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1200 {
1201 gmx_fatal(FARGS,
1202 "The [molecules] section of your topology specifies more than one block of\n"
1203 "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1204 "are trying to partition your solvent into different *groups* (e.g. for\n"
1205 "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1206 "files specify groups. Otherwise, you may wish to change the least-used\n"
1207 "block of molecules with SETTLE constraints into 3 normal constraints.");
1208 }
1209 }
1210 }
1211
1212 constr->n_at2settle_mt = mtop->nmoltype;
1213 snew(constr->at2settle_mt,constr->n_at2settle_mt);
1214 for(mt=0; mt<mtop->nmoltype; mt++)
1215 {
1216 constr->at2settle_mt[mt] =
1217 make_at2settle(mtop->moltype[mt].atoms.nr,
1218 &mtop->moltype[mt].ilist[F_SETTLE]);
1219 }
1220 }
1221
1222 constr->maxwarn = 999;
1223 env = getenv("GMX_MAXCONSTRWARN");
1224 if (env)
1225 {
1226 constr->maxwarn = 0;
1227 sscanf(env,"%d",&constr->maxwarn);
1228 if (fplog)
1229 {
1230 fprintf(fplog,
1231 "Setting the maximum number of constraint warnings to %d\n",
1232 constr->maxwarn);
1233 }
1234 if (MASTER(cr))
1235 {
1236 fprintf(stderr,
1237 "Setting the maximum number of constraint warnings to %d\n",
1238 constr->maxwarn);
1239 }
1240 }
1241 if (constr->maxwarn < 0 && fplog)
1242 {
1243 fprintf(fplog,"maxwarn < 0, will not stop on constraint errors\n");
1244 }
1245 constr->warncount_lincs = 0;
1246 constr->warncount_settle = 0;
1247
1248 /* Initialize the essential dynamics sampling.
1249 * Put the pointer to the ED struct in constr */
1250 constr->ed = ed;
1251 if (ed != NULL)
1252 {
1253 init_edsam(mtop,ir,cr,ed,state->x,state->box);
1254 }
1255
1256 constr->warn_mtop = mtop;
1257
1258 return constr;
1259 }
1260
1261 const t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1262 {
1263 return constr->at2con_mt;
1264 }
1265
1266 const int **atom2settle_moltype(gmx_constr_t constr)
1267 {
1268 return (const int **)constr->at2settle_mt;
1269 }
1270
1271
1272 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
1273 {
1274 const gmx_moltype_t *molt;
1275 const t_block *cgs;
1276 const t_ilist *il;
1277 int mb;
1278 int nat,*at2cg,cg,a,ftype,i;
1279 gmx_bool bInterCG;
1280
1281 bInterCG = FALSE;
1282 for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1283 {
1284 molt = &mtop->moltype[mtop->molblock[mb].type];
1285
1286 if (molt->ilist[F_CONSTR].nr > 0 ||
1287 molt->ilist[F_CONSTRNC].nr > 0 ||
1288 molt->ilist[F_SETTLE].nr > 0)
1289 {
1290 cgs = &molt->cgs;
1291 snew(at2cg,molt->atoms.nr);
1292 for(cg=0; cg<cgs->nr; cg++)
1293 {
1294 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1295 at2cg[a] = cg;
1296 }
1297
1298 for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++)
1299 {
1300 il = &molt->ilist[ftype];
1301 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(ftype))
1302 {
1303 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1304 {
1305 bInterCG = TRUE;
1306 }
1307 }
1308 }
1309
1310 sfree(at2cg);
1311 }
1312 }
1313
1314 return bInterCG;
1315 }
1316
1317 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop)
1318 {
1319 const gmx_moltype_t *molt;
1320 const t_block *cgs;
1321 const t_ilist *il;
1322 int mb;
1323 int nat,*at2cg,cg,a,ftype,i;
1324 gmx_bool bInterCG;
1325
1326 bInterCG = FALSE;
1327 for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1328 {
1329 molt = &mtop->moltype[mtop->molblock[mb].type];
1330
1331 if (molt->ilist[F_SETTLE].nr > 0)
1332 {
1333 cgs = &molt->cgs;
1334 snew(at2cg,molt->atoms.nr);
1335 for(cg=0; cg<cgs->nr; cg++)
1336 {
1337 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1338 at2cg[a] = cg;
1339 }
1340
1341 for(ftype=F_SETTLE; ftype<=F_SETTLE; ftype++)
1342 {
1343 il = &molt->ilist[ftype];
1344 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(F_SETTLE))
1345 {
1346 if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1347 at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1348 {
1349 bInterCG = TRUE;
1350 }
1351 }
1352 }
1353
1354 sfree(at2cg);
1355 }
1356 }
1357
1358 return bInterCG;
1359 }
1360
1361 /* helper functions for andersen temperature control, because the
1362 * gmx_constr construct is only defined in constr.c. Return the list
1363 * of blocks (get_sblock) and the number of blocks (get_nblocks). */
1364
1365 extern int *get_sblock(struct gmx_constr *constr)
1366 {
1367 return constr->sblock;
1368 }
1369
1370 extern int get_nblocks(struct gmx_constr *constr)
1371 {
1372 return constr->nblocks;
1373 }
The ultimate molecular dynamics simulation package