1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
36 /* This file is completely threadsafe - keep it that way! */
53 #include "mtop_util.h"
56 typedef struct gmx_lincsdata {
57 int ncg; /* the global number of constraints */
58 int ncg_flex; /* the global number of flexible constraints */
59 int ncg_triangle;/* the global number of constraints in triangles */
60 int nIter; /* the number of iterations */
61 int nOrder; /* the order of the matrix expansion */
62 int nc; /* the number of constraints */
63 int nc_alloc; /* the number we allocated memory for */
64 int ncc; /* the number of constraint connections */
65 int ncc_alloc; /* the number we allocated memory for */
66 real matlam; /* the FE lambda value used for filling blc and blmf */
67 real *bllen0; /* the reference distance in topology A */
68 real *ddist; /* the reference distance in top B - the r.d. in top A */
69 int *bla; /* the atom pairs involved in the constraints */
70 real *blc; /* 1/sqrt(invmass1 + invmass2) */
71 real *blc1; /* as blc, but with all masses 1 */
72 int *blnr; /* index into blbnb and blmf */
73 int *blbnb; /* list of constraint connections */
74 int ntriangle; /* the local number of constraints in triangles */
75 int *triangle; /* the list of triangle constraints */
76 int *tri_bits; /* the bits tell if the matrix element should be used */
77 int ncc_triangle;/* the number of constraint connections in triangles */
78 real *blmf; /* matrix of mass factors for constraint connections */
79 real *blmf1; /* as blmf, but with all masses 1 */
80 real *bllen; /* the reference bond length */
81 /* arrays for temporary storage in the LINCS algorithm */
87 real *lambda; /* the Lagrange multipliers */
88 /* storage for the constraint RMS relative deviation output */
92 real *lincs_rmsd_data(struct gmx_lincsdata *lincsd)
94 return lincsd->rmsd_data;
97 real lincs_rmsd(struct gmx_lincsdata *lincsd,gmx_bool bSD2)
99 if (lincsd->rmsd_data[0] > 0)
101 return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]);
109 static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd,
111 real *rhs1,real *rhs2,real *sol)
113 int nrec,rec,ncons,b,j,n,nr0,nr1;
115 int ntriangle,tb,bits;
116 const int *blnr=lincsd->blnr,*blbnb=lincsd->blbnb;
117 const int *triangle=lincsd->triangle,*tri_bits=lincsd->tri_bits;
120 ntriangle = lincsd->ntriangle;
121 nrec = lincsd->nOrder;
123 for(rec=0; rec<nrec; rec++)
125 for(b=0; b<ncons; b++)
128 for(n=blnr[b]; n<blnr[b+1]; n++)
131 mvb = mvb + blcc[n]*rhs1[j];
134 sol[b] = sol[b] + mvb;
139 } /* nrec*(ncons+2*nrtot) flops */
143 /* Perform an extra nrec recursions for only the constraints
144 * involved in rigid triangles.
145 * In this way their accuracy should come close to those of the other
146 * constraints, since traingles of constraints can produce eigenvalues
147 * around 0.7, while the effective eigenvalue for bond constraints
148 * is around 0.4 (and 0.7*0.7=0.5).
150 /* We need to copy the temporary array, since only the elements
151 * for constraints involved in triangles are updated
152 * and then the pointers are swapped.
154 for(b=0; b<ncons; b++)
158 for(rec=0; rec<nrec; rec++)
160 for(tb=0; tb<ntriangle; tb++)
167 for(n=nr0; n<nr1; n++)
169 if (bits & (1<<(n-nr0)))
172 mvb = mvb + blcc[n]*rhs1[j];
176 sol[b] = sol[b] + mvb;
181 } /* flops count is missing here */
185 static void do_lincsp(rvec *x,rvec *f,rvec *fp,t_pbc *pbc,
186 struct gmx_lincsdata *lincsd,real *invmass,
187 int econq,real *dvdlambda,
188 gmx_bool bCalcVir,tensor rmdf)
191 real tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,wfac,lam;
193 int ncons,*bla,*blnr,*blbnb;
195 real *blc,*blmf,*blcc,*rhs1,*rhs2,*sol;
201 blbnb = lincsd->blbnb;
202 if (econq != econqForce)
204 /* Use mass-weighted parameters */
210 /* Use non mass-weighted parameters */
212 blmf = lincsd->blmf1;
214 blcc = lincsd->tmpncc;
219 if (econq != econqForce)
224 /* Compute normalized i-j vectors */
227 for(b=0; b<ncons; b++)
229 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
235 for(b=0; b<ncons; b++)
237 rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
239 } /* 16 ncons flops */
242 for(b=0; b<ncons; b++)
249 for(n=blnr[b]; n<blnr[b+1]; n++)
252 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
254 mvb = blc[b]*(tmp0*(f[i][0] - f[j][0]) +
255 tmp1*(f[i][1] - f[j][1]) +
256 tmp2*(f[i][2] - f[j][2]));
261 /* Together: 23*ncons + 6*nrtot flops */
263 lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
264 /* nrec*(ncons+2*nrtot) flops */
266 if (econq != econqForce)
268 for(b=0; b<ncons; b++)
270 /* With econqDeriv_FlexCon only use the flexible constraints */
271 if (econq != econqDeriv_FlexCon ||
272 (lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0))
282 fp[i][0] -= tmp0*im1;
283 fp[i][1] -= tmp1*im1;
284 fp[i][2] -= tmp2*im1;
285 fp[j][0] += tmp0*im2;
286 fp[j][1] += tmp1*im2;
287 fp[j][2] += tmp2*im2;
290 /* This is only correct with forces and invmass=1 */
291 *dvdlambda -= mvb*lincsd->ddist[b];
294 } /* 16 ncons flops */
298 for(b=0; b<ncons; b++)
314 *dvdlambda -= mvb*lincsd->ddist[b];
322 /* Constraint virial,
323 * determines sum r_bond x delta f,
324 * where delta f is the constraint correction
325 * of the quantity that is being constrained.
327 for(b=0; b<ncons; b++)
329 mvb = lincsd->bllen[b]*blc[b]*sol[b];
335 rmdf[i][j] += tmp1*r[b][j];
338 } /* 23 ncons flops */
342 static void do_lincs(rvec *x,rvec *xp,matrix box,t_pbc *pbc,
343 struct gmx_lincsdata *lincsd,real *invmass,
345 real wangle,int *warn,
347 gmx_bool bCalcVir,tensor rmdr)
350 real tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,len2,dlen2,wfac,lam;
352 int ncons,*bla,*blnr,*blbnb;
354 real *blc,*blmf,*bllen,*blcc,*rhs1,*rhs2,*sol,*lambda;
361 blbnb = lincsd->blbnb;
364 bllen = lincsd->bllen;
365 blcc = lincsd->tmpncc;
369 lambda = lincsd->lambda;
371 if (DOMAINDECOMP(cr) && cr->dd->constraints)
373 nlocat = dd_constraints_nlocalatoms(cr->dd);
375 else if (PARTDECOMP(cr))
377 nlocat = pd_constraints_nlocalatoms(cr->pd);
388 /* Compute normalized i-j vectors */
389 for(b=0; b<ncons; b++)
391 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
394 for(b=0; b<ncons; b++)
396 for(n=blnr[b]; n<blnr[b+1]; n++)
398 blcc[n] = blmf[n]*iprod(r[b],r[blbnb[n]]);
400 pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
401 mvb = blc[b]*(iprod(r[b],dx) - bllen[b]);
408 /* Compute normalized i-j vectors */
409 for(b=0; b<ncons; b++)
413 tmp0 = x[i][0] - x[j][0];
414 tmp1 = x[i][1] - x[j][1];
415 tmp2 = x[i][2] - x[j][2];
416 rlen = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
420 } /* 16 ncons flops */
422 for(b=0; b<ncons; b++)
430 for(n=blnr[b]; n<blnr[b+1]; n++)
433 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
435 mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) +
436 tmp1*(xp[i][1] - xp[j][1]) +
437 tmp2*(xp[i][2] - xp[j][2]) - len);
442 /* Together: 26*ncons + 6*nrtot flops */
445 lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
446 /* nrec*(ncons+2*nrtot) flops */
448 for(b=0; b<ncons; b++)
459 xp[i][0] -= tmp0*im1;
460 xp[i][1] -= tmp1*im1;
461 xp[i][2] -= tmp2*im1;
462 xp[j][0] += tmp0*im2;
463 xp[j][1] += tmp1*im2;
464 xp[j][2] += tmp2*im2;
465 } /* 16 ncons flops */
469 ******** Correction for centripetal effects ********
472 wfac = cos(DEG2RAD*wangle);
475 for(iter=0; iter<lincsd->nIter; iter++)
477 if (DOMAINDECOMP(cr) && cr->dd->constraints)
479 /* Communicate the corrected non-local coordinates */
480 dd_move_x_constraints(cr->dd,box,xp,NULL);
482 else if (PARTDECOMP(cr))
484 pd_move_x_constraints(cr,xp,NULL);
487 for(b=0; b<ncons; b++)
492 pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
496 rvec_sub(xp[bla[2*b]],xp[bla[2*b+1]],dx);
499 dlen2 = 2*len2 - norm2(dx);
500 if (dlen2 < wfac*len2 && (nlocat==NULL || nlocat[b]))
506 mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2));
514 } /* 20*ncons flops */
516 lincs_matrix_expand(lincsd,blcc,rhs1,rhs2,sol);
517 /* nrec*(ncons+2*nrtot) flops */
519 for(b=0; b<ncons; b++)
525 lambda[b] = lam - mvb;
531 xp[i][0] -= tmp0*im1;
532 xp[i][1] -= tmp1*im1;
533 xp[i][2] -= tmp2*im1;
534 xp[j][0] += tmp0*im2;
535 xp[j][1] += tmp1*im2;
536 xp[j][2] += tmp2*im2;
537 } /* 17 ncons flops */
538 } /* nit*ncons*(37+9*nrec) flops */
542 /* Correct the velocities */
543 for(b=0; b<ncons; b++)
547 im1 = invmass[i]*lambda[b]*invdt;
548 im2 = invmass[j]*lambda[b]*invdt;
549 v[i][0] += im1*r[b][0];
550 v[i][1] += im1*r[b][1];
551 v[i][2] += im1*r[b][2];
552 v[j][0] -= im2*r[b][0];
553 v[j][1] -= im2*r[b][1];
554 v[j][2] -= im2*r[b][2];
555 } /* 16 ncons flops */
560 /* Only account for local atoms */
561 for(b=0; b<ncons; b++)
563 lambda[b] *= 0.5*nlocat[b];
569 /* Constraint virial */
570 for(b=0; b<ncons; b++)
572 tmp0 = bllen[b]*lambda[b];
578 rmdr[i][j] -= tmp1*r[b][j];
581 } /* 22 ncons flops */
585 * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
586 * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
588 * (26+nrec)*ncons + (6+2*nrec)*nrtot
589 * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
591 * (63+nrec)*ncons + (6+4*nrec)*nrtot
595 void set_lincs_matrix(struct gmx_lincsdata *li,real *invmass,real lambda)
597 int i,a1,a2,n,k,sign,center;
599 const real invsqrt2=0.7071067811865475244;
601 for(i=0; (i<li->nc); i++)
605 li->blc[i] = gmx_invsqrt(invmass[a1] + invmass[a2]);
606 li->blc1[i] = invsqrt2;
609 /* Construct the coupling coefficient matrix blmf */
611 li->ncc_triangle = 0;
612 for(i=0; (i<li->nc); i++)
616 for(n=li->blnr[i]; (n<li->blnr[i+1]); n++)
619 if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
627 if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
637 li->blmf[n] = sign*invmass[center]*li->blc[i]*li->blc[k];
638 li->blmf1[n] = sign*0.5;
639 if (li->ncg_triangle > 0)
641 /* Look for constraint triangles */
642 for(nk=li->blnr[k]; (nk<li->blnr[k+1]); nk++)
645 if (kk != i && kk != k &&
646 (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
648 if (li->ntriangle == 0 ||
649 li->triangle[li->ntriangle-1] < i)
651 /* Add this constraint to the triangle list */
652 li->triangle[li->ntriangle] = i;
653 li->tri_bits[li->ntriangle] = 0;
655 if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
657 gmx_fatal(FARGS,"A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
658 li->blnr[i+1] - li->blnr[i],
659 sizeof(li->tri_bits[0])*8-1);
662 li->tri_bits[li->ntriangle-1] |= (1<<(n-li->blnr[i]));
672 fprintf(debug,"Of the %d constraints %d participate in triangles\n",
673 li->nc,li->ntriangle);
674 fprintf(debug,"There are %d couplings of which %d in triangles\n",
675 li->ncc,li->ncc_triangle);
679 * so we know with which lambda value the masses have been set.
684 static int count_triangle_constraints(t_ilist *ilist,t_blocka *at2con)
687 int c0,a00,a01,n1,c1,a10,a11,ac1,n2,c2,a20,a21;
690 t_iatom *ia1,*ia2,*iap;
692 ncon1 = ilist[F_CONSTR].nr/3;
693 ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
695 ia1 = ilist[F_CONSTR].iatoms;
696 ia2 = ilist[F_CONSTRNC].iatoms;
699 for(c0=0; c0<ncon_tot; c0++)
702 iap = constr_iatomptr(ncon1,ia1,ia2,c0);
705 for(n1=at2con->index[a01]; n1<at2con->index[a01+1]; n1++)
710 iap = constr_iatomptr(ncon1,ia1,ia2,c1);
721 for(n2=at2con->index[ac1]; n2<at2con->index[ac1+1]; n2++)
724 if (c2 != c0 && c2 != c1)
726 iap = constr_iatomptr(ncon1,ia1,ia2,c2);
729 if (a20 == a00 || a21 == a00)
743 return ncon_triangle;
746 static int int_comp(const void *a,const void *b)
748 return (*(int *)a) - (*(int *)b);
751 gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
752 int nflexcon_global,t_blocka *at2con,
753 gmx_bool bPLINCS,int nIter,int nProjOrder)
755 struct gmx_lincsdata *li;
761 fprintf(fplog,"\nInitializing%s LINear Constraint Solver\n",
762 bPLINCS ? " Parallel" : "");
768 gmx_mtop_ftype_count(mtop,F_CONSTR) +
769 gmx_mtop_ftype_count(mtop,F_CONSTRNC);
770 li->ncg_flex = nflexcon_global;
772 li->ncg_triangle = 0;
773 for(mb=0; mb<mtop->nmolblock; mb++)
775 molt = &mtop->moltype[mtop->molblock[mb].type];
777 mtop->molblock[mb].nmol*
778 count_triangle_constraints(molt->ilist,
779 &at2con[mtop->molblock[mb].type]);
783 li->nOrder = nProjOrder;
785 if (bPLINCS || li->ncg_triangle > 0)
787 please_cite(fplog,"Hess2008a");
791 please_cite(fplog,"Hess97a");
796 fprintf(fplog,"The number of constraints is %d\n",li->ncg);
799 fprintf(fplog,"There are inter charge-group constraints,\n"
800 "will communicate selected coordinates each lincs iteration\n");
802 if (li->ncg_triangle > 0)
805 "%d constraints are involved in constraint triangles,\n"
806 "will apply an additional matrix expansion of order %d for couplings\n"
807 "between constraints inside triangles\n",
808 li->ncg_triangle,li->nOrder);
815 void set_lincs(t_idef *idef,t_mdatoms *md,
816 gmx_bool bDynamics,t_commrec *cr,
817 struct gmx_lincsdata *li)
819 int start,natoms,nflexcon;
822 int i,k,ncc_alloc,ni,con,nconnect,concon;
830 /* This is the local topology, so there are only F_CONSTR constraints */
831 if (idef->il[F_CONSTR].nr == 0)
833 /* There are no constraints,
834 * we do not need to fill any data structures.
841 fprintf(debug,"Building the LINCS connectivity\n");
844 if (DOMAINDECOMP(cr))
846 if (cr->dd->constraints)
848 dd_get_constraint_range(cr->dd,&start,&natoms);
852 natoms = cr->dd->nat_home;
856 else if(PARTDECOMP(cr))
858 pd_get_constraint_range(cr->pd,&start,&natoms);
865 at2con = make_at2con(start,natoms,idef->il,idef->iparams,bDynamics,
869 if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)
871 li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);
872 srenew(li->bllen0,li->nc_alloc);
873 srenew(li->ddist,li->nc_alloc);
874 srenew(li->bla,2*li->nc_alloc);
875 srenew(li->blc,li->nc_alloc);
876 srenew(li->blc1,li->nc_alloc);
877 srenew(li->blnr,li->nc_alloc+1);
878 srenew(li->bllen,li->nc_alloc);
879 srenew(li->tmpv,li->nc_alloc);
880 srenew(li->tmp1,li->nc_alloc);
881 srenew(li->tmp2,li->nc_alloc);
882 srenew(li->tmp3,li->nc_alloc);
883 srenew(li->lambda,li->nc_alloc);
884 if (li->ncg_triangle > 0)
886 /* This is allocating too much, but it is difficult to improve */
887 srenew(li->triangle,li->nc_alloc);
888 srenew(li->tri_bits,li->nc_alloc);
892 iatom = idef->il[F_CONSTR].iatoms;
894 ncc_alloc = li->ncc_alloc;
897 ni = idef->il[F_CONSTR].nr/3;
901 li->blnr[con] = nconnect;
908 lenA = idef->iparams[type].constr.dA;
909 lenB = idef->iparams[type].constr.dB;
910 /* Skip the flexible constraints when not doing dynamics */
911 if (bDynamics || lenA!=0 || lenB!=0)
913 li->bllen0[con] = lenA;
914 li->ddist[con] = lenB - lenA;
915 /* Set the length to the topology A length */
916 li->bllen[con] = li->bllen0[con];
918 li->bla[2*con+1] = a2;
919 /* Construct the constraint connection matrix blbnb */
920 for(k=at2con.index[a1-start]; k<at2con.index[a1-start+1]; k++)
922 concon = at2con.a[k];
925 if (nconnect >= ncc_alloc)
927 ncc_alloc = over_alloc_small(nconnect+1);
928 srenew(li->blbnb,ncc_alloc);
930 li->blbnb[nconnect++] = concon;
933 for(k=at2con.index[a2-start]; k<at2con.index[a2-start+1]; k++)
935 concon = at2con.a[k];
938 if (nconnect+1 > ncc_alloc)
940 ncc_alloc = over_alloc_small(nconnect+1);
941 srenew(li->blbnb,ncc_alloc);
943 li->blbnb[nconnect++] = concon;
946 li->blnr[con+1] = nconnect;
950 /* Order the blbnb matrix to optimize memory access */
951 qsort(&(li->blbnb[li->blnr[con]]),li->blnr[con+1]-li->blnr[con],
952 sizeof(li->blbnb[0]),int_comp);
954 /* Increase the constraint count */
959 done_blocka(&at2con);
961 /* This is the real number of constraints,
962 * without dynamics the flexible constraints are not present.
966 li->ncc = li->blnr[con];
969 /* Since the matrix is static, we can free some memory */
971 srenew(li->blbnb,ncc_alloc);
974 if (ncc_alloc > li->ncc_alloc)
976 li->ncc_alloc = ncc_alloc;
977 srenew(li->blmf,li->ncc_alloc);
978 srenew(li->blmf1,li->ncc_alloc);
979 srenew(li->tmpncc,li->ncc_alloc);
984 fprintf(debug,"Number of constraints is %d, couplings %d\n",
988 set_lincs_matrix(li,md->invmass,md->lambda);
991 static void lincs_warning(FILE *fplog,
992 gmx_domdec_t *dd,rvec *x,rvec *xprime,t_pbc *pbc,
993 int ncons,int *bla,real *bllen,real wangle,
994 int maxwarn,int *warncount)
998 real wfac,d0,d1,cosine;
1001 wfac=cos(DEG2RAD*wangle);
1003 sprintf(buf,"bonds that rotated more than %g degrees:\n"
1004 " atom 1 atom 2 angle previous, current, constraint length\n",
1006 fprintf(stderr,"%s",buf);
1009 fprintf(fplog,"%s",buf);
1012 for(b=0;b<ncons;b++)
1018 pbc_dx_aiuc(pbc,x[i],x[j],v0);
1019 pbc_dx_aiuc(pbc,xprime[i],xprime[j],v1);
1023 rvec_sub(x[i],x[j],v0);
1024 rvec_sub(xprime[i],xprime[j],v1);
1028 cosine = iprod(v0,v1)/(d0*d1);
1031 sprintf(buf," %6d %6d %5.1f %8.4f %8.4f %8.4f\n",
1032 ddglatnr(dd,i),ddglatnr(dd,j),
1033 RAD2DEG*acos(cosine),d0,d1,bllen[b]);
1034 fprintf(stderr,"%s",buf);
1037 fprintf(fplog,"%s",buf);
1039 if (!gmx_isfinite(d1))
1041 gmx_fatal(FARGS,"Bond length not finite.");
1047 if (*warncount > maxwarn)
1049 too_many_constraint_warnings(econtLINCS,*warncount);
1053 static void cconerr(gmx_domdec_t *dd,
1054 int ncons,int *bla,real *bllen,rvec *x,t_pbc *pbc,
1055 real *ncons_loc,real *ssd,real *max,int *imax)
1057 real len,d,ma,ssd2,r2;
1058 int *nlocat,count,b,im;
1061 if (dd && dd->constraints)
1063 nlocat = dd_constraints_nlocalatoms(dd);
1074 for(b=0;b<ncons;b++)
1078 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
1081 rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
1084 len = r2*gmx_invsqrt(r2);
1085 d = fabs(len/bllen[b]-1);
1086 if (d > ma && (nlocat==NULL || nlocat[b]))
1098 ssd2 += nlocat[b]*d*d;
1103 *ncons_loc = (nlocat ? 0.5 : 1)*count;
1104 *ssd = (nlocat ? 0.5 : 1)*ssd2;
1109 static void dump_conf(gmx_domdec_t *dd,struct gmx_lincsdata *li,
1111 char *name,gmx_bool bAll,rvec *x,matrix box)
1117 dd_get_constraint_range(dd,&ac0,&ac1);
1119 sprintf(str,"%s_%d_%d_%d.pdb",name,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
1120 fp = gmx_fio_fopen(str,"w");
1121 fprintf(fp,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1\n",
1122 10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
1124 for(i=0; i<ac1; i++)
1126 if (i < dd->nat_home || (bAll && i >= ac0 && i < ac1))
1128 fprintf(fp,"%-6s%5u %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
1129 "ATOM",ddglatnr(dd,i),"C","ALA",' ',i+1,
1130 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],
1131 1.0,i<dd->nat_tot ? 0.0 : 1.0);
1136 for(i=0; i<li->nc; i++)
1138 fprintf(fp,"CONECT%5d%5d\n",
1139 ddglatnr(dd,li->bla[2*i]),
1140 ddglatnr(dd,li->bla[2*i+1]));
1146 gmx_bool constrain_lincs(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
1148 gmx_large_int_t step,
1149 struct gmx_lincsdata *lincsd,t_mdatoms *md,
1151 rvec *x,rvec *xprime,rvec *min_proj,matrix box,
1152 real lambda,real *dvdlambda,
1154 gmx_bool bCalcVir,tensor rmdr,
1157 int maxwarn,int *warncount)
1159 char buf[STRLEN],buf2[22],buf3[STRLEN];
1160 int i,warn,p_imax,error;
1161 real ncons_loc,p_ssd,p_max=0;
1162 t_pbc pbc,*pbc_null;
1168 if (lincsd->nc == 0 && cr->dd == NULL)
1172 lincsd->rmsd_data[0] = 0;
1173 if (ir->eI == eiSD2 && v == NULL)
1181 lincsd->rmsd_data[i] = 0;
1187 /* We do not need full pbc when constraints do not cross charge groups,
1188 * i.e. when dd->constraint_comm==NULL
1190 if ((cr->dd || ir->bPeriodicMols) && !(cr->dd && cr->dd->constraint_comm==NULL))
1192 /* With pbc=screw the screw has been changed to a shift
1193 * by the constraint coordinate communication routine,
1194 * so that here we can use normal pbc.
1196 pbc_null = set_pbc_dd(&pbc,ir->ePBC,cr->dd,FALSE,box);
1204 /* Communicate the coordinates required for the non-local constraints */
1205 dd_move_x_constraints(cr->dd,box,x,xprime);
1206 /* dump_conf(dd,lincsd,NULL,"con",TRUE,xprime,box); */
1208 else if (PARTDECOMP(cr))
1210 pd_move_x_constraints(cr,x,xprime);
1213 if (econq == econqCoord)
1215 if (ir->efep != efepNO)
1217 if (md->nMassPerturbed && lincsd->matlam != md->lambda)
1219 set_lincs_matrix(lincsd,md->invmass,md->lambda);
1222 for(i=0; i<lincsd->nc; i++)
1224 lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
1228 if (lincsd->ncg_flex)
1230 /* Set the flexible constraint lengths to the old lengths */
1233 for(i=0; i<lincsd->nc; i++)
1235 if (lincsd->bllen[i] == 0) {
1236 pbc_dx_aiuc(pbc_null,x[lincsd->bla[2*i]],x[lincsd->bla[2*i+1]],dx);
1237 lincsd->bllen[i] = norm(dx);
1243 for(i=0; i<lincsd->nc; i++)
1245 if (lincsd->bllen[i] == 0)
1248 sqrt(distance2(x[lincsd->bla[2*i]],
1249 x[lincsd->bla[2*i+1]]));
1257 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1258 &ncons_loc,&p_ssd,&p_max,&p_imax);
1261 do_lincs(x,xprime,box,pbc_null,lincsd,md->invmass,cr,
1262 ir->LincsWarnAngle,&warn,
1263 invdt,v,bCalcVir,rmdr);
1265 if (ir->efep != efepNO)
1269 dt_2 = 1.0/(ir->delta_t*ir->delta_t);
1270 for(i=0; (i<lincsd->nc); i++)
1272 dvdl += lincsd->lambda[i]*dt_2*lincsd->ddist[i];
1277 if (bLog && fplog && lincsd->nc > 0)
1279 fprintf(fplog," Rel. Constraint Deviation: RMS MAX between atoms\n");
1280 fprintf(fplog," Before LINCS %.6f %.6f %6d %6d\n",
1281 sqrt(p_ssd/ncons_loc),p_max,
1282 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1283 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1287 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1288 &ncons_loc,&p_ssd,&p_max,&p_imax);
1289 /* Check if we are doing the second part of SD */
1290 if (ir->eI == eiSD2 && v == NULL)
1298 lincsd->rmsd_data[0] = ncons_loc;
1299 lincsd->rmsd_data[i] = p_ssd;
1303 lincsd->rmsd_data[0] = 0;
1304 lincsd->rmsd_data[1] = 0;
1305 lincsd->rmsd_data[2] = 0;
1307 if (bLog && fplog && lincsd->nc > 0)
1310 " After LINCS %.6f %.6f %6d %6d\n\n",
1311 sqrt(p_ssd/ncons_loc),p_max,
1312 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1313 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1320 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc_null,
1321 &ncons_loc,&p_ssd,&p_max,&p_imax);
1324 sprintf(buf3," in simulation %d", cr->ms->sim);
1330 sprintf(buf,"\nStep %s, time %g (ps) LINCS WARNING%s\n"
1331 "relative constraint deviation after LINCS:\n"
1332 "rms %.6f, max %.6f (between atoms %d and %d)\n",
1333 gmx_step_str(step,buf2),ir->init_t+step*ir->delta_t,
1335 sqrt(p_ssd/ncons_loc),p_max,
1336 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1337 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1340 fprintf(fplog,"%s",buf);
1342 fprintf(stderr,"%s",buf);
1343 lincs_warning(fplog,cr->dd,x,xprime,pbc_null,
1344 lincsd->nc,lincsd->bla,lincsd->bllen,
1345 ir->LincsWarnAngle,maxwarn,warncount);
1347 bOK = (p_max < 0.5);
1350 if (lincsd->ncg_flex) {
1351 for(i=0; (i<lincsd->nc); i++)
1352 if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
1353 lincsd->bllen[i] = 0;
1358 do_lincsp(x,xprime,min_proj,pbc_null,lincsd,md->invmass,econq,dvdlambda,
1362 /* count assuming nit=1 */
1363 inc_nrnb(nrnb,eNR_LINCS,lincsd->nc);
1364 inc_nrnb(nrnb,eNR_LINCSMAT,(2+lincsd->nOrder)*lincsd->ncc);
1365 if (lincsd->ntriangle > 0)
1367 inc_nrnb(nrnb,eNR_LINCSMAT,lincsd->nOrder*lincsd->ncc_triangle);
1371 inc_nrnb(nrnb,eNR_CONSTR_V,lincsd->nc*2);
1375 inc_nrnb(nrnb,eNR_CONSTR_VIR,lincsd->nc);