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"
55 #include "gmx_omp_nthreads.h"
59 int b0; /* first constraint for this thread */
60 int b1; /* b1-1 is the last constraint for this thread */
61 int nind; /* number of indices */
62 int *ind; /* constraint index for updating atom data */
63 int nind_r; /* number of indices */
64 int *ind_r; /* constraint index for updating atom data */
65 int ind_nalloc; /* allocation size of ind and ind_r */
66 tensor vir_r_m_dr;/* temporary variable for virial calculation */
69 typedef struct gmx_lincsdata {
70 int ncg; /* the global number of constraints */
71 int ncg_flex; /* the global number of flexible constraints */
72 int ncg_triangle;/* the global number of constraints in triangles */
73 int nIter; /* the number of iterations */
74 int nOrder; /* the order of the matrix expansion */
75 int nc; /* the number of constraints */
76 int nc_alloc; /* the number we allocated memory for */
77 int ncc; /* the number of constraint connections */
78 int ncc_alloc; /* the number we allocated memory for */
79 real matlam; /* the FE lambda value used for filling blc and blmf */
80 real *bllen0; /* the reference distance in topology A */
81 real *ddist; /* the reference distance in top B - the r.d. in top A */
82 int *bla; /* the atom pairs involved in the constraints */
83 real *blc; /* 1/sqrt(invmass1 + invmass2) */
84 real *blc1; /* as blc, but with all masses 1 */
85 int *blnr; /* index into blbnb and blmf */
86 int *blbnb; /* list of constraint connections */
87 int ntriangle; /* the local number of constraints in triangles */
88 int *triangle; /* the list of triangle constraints */
89 int *tri_bits; /* the bits tell if the matrix element should be used */
90 int ncc_triangle;/* the number of constraint connections in triangles */
91 real *blmf; /* matrix of mass factors for constraint connections */
92 real *blmf1; /* as blmf, but with all masses 1 */
93 real *bllen; /* the reference bond length */
94 int nth; /* The number of threads doing LINCS */
95 lincs_thread_t *th; /* LINCS thread division */
96 unsigned *atf; /* atom flags for thread parallelization */
97 int atf_nalloc; /* allocation size of atf */
98 /* arrays for temporary storage in the LINCS algorithm */
105 real *mlambda; /* the Lagrange multipliers * -1 */
106 /* storage for the constraint RMS relative deviation output */
110 real *lincs_rmsd_data(struct gmx_lincsdata *lincsd)
112 return lincsd->rmsd_data;
115 real lincs_rmsd(struct gmx_lincsdata *lincsd,gmx_bool bSD2)
117 if (lincsd->rmsd_data[0] > 0)
119 return sqrt(lincsd->rmsd_data[bSD2 ? 2 : 1]/lincsd->rmsd_data[0]);
127 /* Do a set of nrec LINCS matrix multiplications.
128 * This function will return with up to date thread-local
129 * constraint data, without an OpenMP barrier.
131 static void lincs_matrix_expand(const struct gmx_lincsdata *lincsd,
134 real *rhs1,real *rhs2,real *sol)
136 int nrec,rec,b,j,n,nr0,nr1;
138 int ntriangle,tb,bits;
139 const int *blnr=lincsd->blnr,*blbnb=lincsd->blbnb;
140 const int *triangle=lincsd->triangle,*tri_bits=lincsd->tri_bits;
142 ntriangle = lincsd->ntriangle;
143 nrec = lincsd->nOrder;
145 for(rec=0; rec<nrec; rec++)
151 for(n=blnr[b]; n<blnr[b+1]; n++)
154 mvb = mvb + blcc[n]*rhs1[j];
157 sol[b] = sol[b] + mvb;
162 } /* nrec*(ncons+2*nrtot) flops */
166 /* Perform an extra nrec recursions for only the constraints
167 * involved in rigid triangles.
168 * In this way their accuracy should come close to those of the other
169 * constraints, since traingles of constraints can produce eigenvalues
170 * around 0.7, while the effective eigenvalue for bond constraints
171 * is around 0.4 (and 0.7*0.7=0.5).
173 /* We need to copy the temporary array, since only the elements
174 * for constraints involved in triangles are updated and then
175 * the pointers are swapped. This saving copying the whole arrary.
176 * We need barrier as other threads might still be reading from rhs2.
186 for(rec=0; rec<nrec; rec++)
188 for(tb=0; tb<ntriangle; tb++)
195 for(n=nr0; n<nr1; n++)
197 if (bits & (1<<(n-nr0)))
200 mvb = mvb + blcc[n]*rhs1[j];
204 sol[b] = sol[b] + mvb;
210 } /* flops count is missing here */
212 /* We need a barrier here as the calling routine will continue
213 * to operate on the thread-local constraints without barrier.
219 static void lincs_update_atoms_noind(int ncons,const int *bla,
221 const real *fac,rvec *r,
226 real mvb,im1,im2,tmp0,tmp1,tmp2;
228 for(b=0; b<ncons; b++)
244 } /* 16 ncons flops */
247 static void lincs_update_atoms_ind(int ncons,const int *ind,const int *bla,
249 const real *fac,rvec *r,
254 real mvb,im1,im2,tmp0,tmp1,tmp2;
256 for(bi=0; bi<ncons; bi++)
273 } /* 16 ncons flops */
276 static void lincs_update_atoms(struct gmx_lincsdata *li,int th,
278 const real *fac,rvec *r,
284 /* Single thread, we simply update for all constraints */
285 lincs_update_atoms_noind(li->nc,li->bla,prefac,fac,r,invmass,x);
289 /* Update the atom vector components for our thread local
290 * constraints that only access our local atom range.
291 * This can be done without a barrier.
293 lincs_update_atoms_ind(li->th[th].nind,li->th[th].ind,
294 li->bla,prefac,fac,r,invmass,x);
296 if (li->th[li->nth].nind > 0)
298 /* Update the constraints that operate on atoms
299 * in multiple thread atom blocks on the master thread.
304 lincs_update_atoms_ind(li->th[li->nth].nind,
306 li->bla,prefac,fac,r,invmass,x);
312 /* LINCS projection, works on derivatives of the coordinates */
313 static void do_lincsp(rvec *x,rvec *f,rvec *fp,t_pbc *pbc,
314 struct gmx_lincsdata *lincsd,int th,
316 int econq,real *dvdlambda,
317 gmx_bool bCalcVir,tensor rmdf)
320 real tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,wfac,lam;
322 int *bla,*blnr,*blbnb;
324 real *blc,*blmf,*blcc,*rhs1,*rhs2,*sol;
326 b0 = lincsd->th[th].b0;
327 b1 = lincsd->th[th].b1;
332 blbnb = lincsd->blbnb;
333 if (econq != econqForce)
335 /* Use mass-weighted parameters */
341 /* Use non mass-weighted parameters */
343 blmf = lincsd->blmf1;
345 blcc = lincsd->tmpncc;
350 /* Compute normalized i-j vectors */
355 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
363 rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
365 } /* 16 ncons flops */
376 for(n=blnr[b]; n<blnr[b+1]; n++)
379 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
381 mvb = blc[b]*(tmp0*(f[i][0] - f[j][0]) +
382 tmp1*(f[i][1] - f[j][1]) +
383 tmp2*(f[i][2] - f[j][2]));
388 /* Together: 23*ncons + 6*nrtot flops */
390 lincs_matrix_expand(lincsd,b0,b1,blcc,rhs1,rhs2,sol);
391 /* nrec*(ncons+2*nrtot) flops */
393 if (econq == econqDeriv_FlexCon)
395 /* We only want to constraint the flexible constraints,
396 * so we mask out the normal ones by setting sol to 0.
400 if (!(lincsd->bllen0[b] == 0 && lincsd->ddist[b] == 0))
407 if (econq != econqForce)
409 lincs_update_atoms(lincsd,th,1.0,sol,r,invmass,fp);
429 if (dvdlambda != NULL)
434 *dvdlambda -= blc[b]*sol[b]*lincsd->ddist[b];
442 /* Constraint virial,
443 * determines sum r_bond x delta f,
444 * where delta f is the constraint correction
445 * of the quantity that is being constrained.
449 mvb = lincsd->bllen[b]*blc[b]*sol[b];
455 rmdf[i][j] += tmp1*r[b][j];
458 } /* 23 ncons flops */
462 static void do_lincs(rvec *x,rvec *xp,matrix box,t_pbc *pbc,
463 struct gmx_lincsdata *lincsd,int th,
466 gmx_bool bCalcLambda,
467 real wangle,int *warn,
469 gmx_bool bCalcVir,tensor vir_r_m_dr)
471 int b0,b1,b,i,j,k,n,iter;
472 real tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,len2,dlen2,wfac;
474 int *bla,*blnr,*blbnb;
476 real *blc,*blmf,*bllen,*blcc,*rhs1,*rhs2,*sol,*blc_sol,*mlambda;
479 b0 = lincsd->th[th].b0;
480 b1 = lincsd->th[th].b1;
485 blbnb = lincsd->blbnb;
488 bllen = lincsd->bllen;
489 blcc = lincsd->tmpncc;
493 blc_sol= lincsd->tmp4;
494 mlambda= lincsd->mlambda;
496 if (DOMAINDECOMP(cr) && cr->dd->constraints)
498 nlocat = dd_constraints_nlocalatoms(cr->dd);
500 else if (PARTDECOMP(cr))
502 nlocat = pd_constraints_nlocalatoms(cr->pd);
511 /* Compute normalized i-j vectors */
514 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
520 for(n=blnr[b]; n<blnr[b+1]; n++)
522 blcc[n] = blmf[n]*iprod(r[b],r[blbnb[n]]);
524 pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
525 mvb = blc[b]*(iprod(r[b],dx) - bllen[b]);
532 /* Compute normalized i-j vectors */
537 tmp0 = x[i][0] - x[j][0];
538 tmp1 = x[i][1] - x[j][1];
539 tmp2 = x[i][2] - x[j][2];
540 rlen = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
544 } /* 16 ncons flops */
555 for(n=blnr[b]; n<blnr[b+1]; n++)
558 blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);
560 mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) +
561 tmp1*(xp[i][1] - xp[j][1]) +
562 tmp2*(xp[i][2] - xp[j][2]) - len);
567 /* Together: 26*ncons + 6*nrtot flops */
570 lincs_matrix_expand(lincsd,b0,b1,blcc,rhs1,rhs2,sol);
571 /* nrec*(ncons+2*nrtot) flops */
575 mlambda[b] = blc[b]*sol[b];
578 /* Update the coordinates */
579 lincs_update_atoms(lincsd,th,1.0,mlambda,r,invmass,xp);
582 ******** Correction for centripetal effects ********
585 wfac = cos(DEG2RAD*wangle);
588 for(iter=0; iter<lincsd->nIter; iter++)
590 if ((DOMAINDECOMP(cr) && cr->dd->constraints) ||
596 /* Communicate the corrected non-local coordinates */
597 if (DOMAINDECOMP(cr))
599 dd_move_x_constraints(cr->dd,box,xp,NULL);
603 pd_move_x_constraints(cr,xp,NULL);
614 pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);
618 rvec_sub(xp[bla[2*b]],xp[bla[2*b+1]],dx);
621 dlen2 = 2*len2 - norm2(dx);
622 if (dlen2 < wfac*len2 && (nlocat==NULL || nlocat[b]))
628 mvb = blc[b]*(len - dlen2*gmx_invsqrt(dlen2));
636 } /* 20*ncons flops */
638 lincs_matrix_expand(lincsd,b0,b1,blcc,rhs1,rhs2,sol);
639 /* nrec*(ncons+2*nrtot) flops */
648 /* Update the coordinates */
649 lincs_update_atoms(lincsd,th,1.0,blc_sol,r,invmass,xp);
651 /* nit*ncons*(37+9*nrec) flops */
655 /* Update the velocities */
656 lincs_update_atoms(lincsd,th,invdt,mlambda,r,invmass,v);
660 if (nlocat != NULL && bCalcLambda)
662 /* In lincs_update_atoms thread might cross-read mlambda */
665 /* Only account for local atoms */
668 mlambda[b] *= 0.5*nlocat[b];
674 /* Constraint virial */
677 tmp0 = -bllen[b]*mlambda[b];
683 vir_r_m_dr[i][j] -= tmp1*r[b][j];
686 } /* 22 ncons flops */
690 * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
691 * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
693 * (26+nrec)*ncons + (6+2*nrec)*nrtot
694 * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
696 * (63+nrec)*ncons + (6+4*nrec)*nrtot
700 void set_lincs_matrix(struct gmx_lincsdata *li,real *invmass,real lambda)
702 int i,a1,a2,n,k,sign,center;
704 const real invsqrt2=0.7071067811865475244;
706 for(i=0; (i<li->nc); i++)
710 li->blc[i] = gmx_invsqrt(invmass[a1] + invmass[a2]);
711 li->blc1[i] = invsqrt2;
714 /* Construct the coupling coefficient matrix blmf */
716 li->ncc_triangle = 0;
717 for(i=0; (i<li->nc); i++)
721 for(n=li->blnr[i]; (n<li->blnr[i+1]); n++)
724 if (a1 == li->bla[2*k] || a2 == li->bla[2*k+1])
732 if (a1 == li->bla[2*k] || a1 == li->bla[2*k+1])
742 li->blmf[n] = sign*invmass[center]*li->blc[i]*li->blc[k];
743 li->blmf1[n] = sign*0.5;
744 if (li->ncg_triangle > 0)
746 /* Look for constraint triangles */
747 for(nk=li->blnr[k]; (nk<li->blnr[k+1]); nk++)
750 if (kk != i && kk != k &&
751 (li->bla[2*kk] == end || li->bla[2*kk+1] == end))
753 if (li->ntriangle == 0 ||
754 li->triangle[li->ntriangle-1] < i)
756 /* Add this constraint to the triangle list */
757 li->triangle[li->ntriangle] = i;
758 li->tri_bits[li->ntriangle] = 0;
760 if (li->blnr[i+1] - li->blnr[i] > sizeof(li->tri_bits[0])*8 - 1)
762 gmx_fatal(FARGS,"A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
763 li->blnr[i+1] - li->blnr[i],
764 sizeof(li->tri_bits[0])*8-1);
767 li->tri_bits[li->ntriangle-1] |= (1<<(n-li->blnr[i]));
777 fprintf(debug,"Of the %d constraints %d participate in triangles\n",
778 li->nc,li->ntriangle);
779 fprintf(debug,"There are %d couplings of which %d in triangles\n",
780 li->ncc,li->ncc_triangle);
784 * so we know with which lambda value the masses have been set.
789 static int count_triangle_constraints(t_ilist *ilist,t_blocka *at2con)
792 int c0,a00,a01,n1,c1,a10,a11,ac1,n2,c2,a20,a21;
795 t_iatom *ia1,*ia2,*iap;
797 ncon1 = ilist[F_CONSTR].nr/3;
798 ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
800 ia1 = ilist[F_CONSTR].iatoms;
801 ia2 = ilist[F_CONSTRNC].iatoms;
804 for(c0=0; c0<ncon_tot; c0++)
807 iap = constr_iatomptr(ncon1,ia1,ia2,c0);
810 for(n1=at2con->index[a01]; n1<at2con->index[a01+1]; n1++)
815 iap = constr_iatomptr(ncon1,ia1,ia2,c1);
826 for(n2=at2con->index[ac1]; n2<at2con->index[ac1+1]; n2++)
829 if (c2 != c0 && c2 != c1)
831 iap = constr_iatomptr(ncon1,ia1,ia2,c2);
834 if (a20 == a00 || a21 == a00)
848 return ncon_triangle;
851 static int int_comp(const void *a,const void *b)
853 return (*(int *)a) - (*(int *)b);
856 gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
857 int nflexcon_global,t_blocka *at2con,
858 gmx_bool bPLINCS,int nIter,int nProjOrder)
860 struct gmx_lincsdata *li;
866 fprintf(fplog,"\nInitializing%s LINear Constraint Solver\n",
867 bPLINCS ? " Parallel" : "");
873 gmx_mtop_ftype_count(mtop,F_CONSTR) +
874 gmx_mtop_ftype_count(mtop,F_CONSTRNC);
875 li->ncg_flex = nflexcon_global;
877 li->ncg_triangle = 0;
878 for(mb=0; mb<mtop->nmolblock; mb++)
880 molt = &mtop->moltype[mtop->molblock[mb].type];
882 mtop->molblock[mb].nmol*
883 count_triangle_constraints(molt->ilist,
884 &at2con[mtop->molblock[mb].type]);
888 li->nOrder = nProjOrder;
890 /* LINCS can run on any number of threads.
891 * Currently the number is fixed for the whole simulation,
892 * but it could be set in set_lincs().
894 li->nth = gmx_omp_nthreads_get(emntLINCS);
901 /* Allocate an extra elements for "thread-overlap" constraints */
902 snew(li->th,li->nth+1);
906 fprintf(debug,"LINCS: using %d threads\n",li->nth);
909 if (bPLINCS || li->ncg_triangle > 0)
911 please_cite(fplog,"Hess2008a");
915 please_cite(fplog,"Hess97a");
920 fprintf(fplog,"The number of constraints is %d\n",li->ncg);
923 fprintf(fplog,"There are inter charge-group constraints,\n"
924 "will communicate selected coordinates each lincs iteration\n");
926 if (li->ncg_triangle > 0)
929 "%d constraints are involved in constraint triangles,\n"
930 "will apply an additional matrix expansion of order %d for couplings\n"
931 "between constraints inside triangles\n",
932 li->ncg_triangle,li->nOrder);
939 /* Sets up the work division over the threads */
940 static void lincs_thread_setup(struct gmx_lincsdata *li,int natoms)
942 lincs_thread_t *li_m;
947 if (natoms > li->atf_nalloc)
949 li->atf_nalloc = over_alloc_large(natoms);
950 srenew(li->atf,li->atf_nalloc);
954 /* Clear the atom flags */
955 for(a=0; a<natoms; a++)
960 for(th=0; th<li->nth; th++)
962 lincs_thread_t *li_th;
967 /* The constraints are divided equally over the threads */
968 li_th->b0 = (li->nc* th )/li->nth;
969 li_th->b1 = (li->nc*(th+1))/li->nth;
971 if (th < sizeof(*atf)*8)
973 /* For each atom set a flag for constraints from each */
974 for(b=li_th->b0; b<li_th->b1; b++)
976 atf[li->bla[b*2] ] |= (1U<<th);
977 atf[li->bla[b*2+1]] |= (1U<<th);
982 #pragma omp parallel for num_threads(li->nth) schedule(static)
983 for(th=0; th<li->nth; th++)
985 lincs_thread_t *li_th;
991 if (li_th->b1 - li_th->b0 > li_th->ind_nalloc)
993 li_th->ind_nalloc = over_alloc_large(li_th->b1-li_th->b0);
994 srenew(li_th->ind,li_th->ind_nalloc);
995 srenew(li_th->ind_r,li_th->ind_nalloc);
998 if (th < sizeof(*atf)*8)
1000 mask = (1U<<th) - 1U;
1004 for(b=li_th->b0; b<li_th->b1; b++)
1006 /* We let the constraint with the lowest thread index
1007 * operate on atoms with constraints from multiple threads.
1009 if (((atf[li->bla[b*2]] & mask) == 0) &&
1010 ((atf[li->bla[b*2+1]] & mask) == 0))
1012 /* Add the constraint to the local atom update index */
1013 li_th->ind[li_th->nind++] = b;
1017 /* Add the constraint to the rest block */
1018 li_th->ind_r[li_th->nind_r++] = b;
1024 /* We are out of bits, assign all constraints to rest */
1025 for(b=li_th->b0; b<li_th->b1; b++)
1027 li_th->ind_r[li_th->nind_r++] = b;
1032 /* We need to copy all constraints which have not be assigned
1033 * to a thread to a separate list which will be handled by one thread.
1035 li_m = &li->th[li->nth];
1038 for(th=0; th<li->nth; th++)
1040 lincs_thread_t *li_th;
1043 li_th = &li->th[th];
1045 if (li_m->nind + li_th->nind_r > li_m->ind_nalloc)
1047 li_m->ind_nalloc = over_alloc_large(li_m->nind+li_th->nind_r);
1048 srenew(li_m->ind,li_m->ind_nalloc);
1051 for(b=0; b<li_th->nind_r; b++)
1053 li_m->ind[li_m->nind++] = li_th->ind_r[b];
1058 fprintf(debug,"LINCS thread %d: %d constraints\n",
1065 fprintf(debug,"LINCS thread r: %d constraints\n",
1071 void set_lincs(t_idef *idef,t_mdatoms *md,
1072 gmx_bool bDynamics,t_commrec *cr,
1073 struct gmx_lincsdata *li)
1075 int start,natoms,nflexcon;
1078 int i,k,ncc_alloc,ni,con,nconnect,concon;
1085 /* Zero the thread index ranges.
1086 * Otherwise without local constraints we could return with old ranges.
1088 for(i=0; i<li->nth; i++)
1096 li->th[li->nth].nind = 0;
1099 /* This is the local topology, so there are only F_CONSTR constraints */
1100 if (idef->il[F_CONSTR].nr == 0)
1102 /* There are no constraints,
1103 * we do not need to fill any data structures.
1110 fprintf(debug,"Building the LINCS connectivity\n");
1113 if (DOMAINDECOMP(cr))
1115 if (cr->dd->constraints)
1117 dd_get_constraint_range(cr->dd,&start,&natoms);
1121 natoms = cr->dd->nat_home;
1125 else if(PARTDECOMP(cr))
1127 pd_get_constraint_range(cr->pd,&start,&natoms);
1132 natoms = md->homenr;
1134 at2con = make_at2con(start,natoms,idef->il,idef->iparams,bDynamics,
1138 if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)
1140 li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);
1141 srenew(li->bllen0,li->nc_alloc);
1142 srenew(li->ddist,li->nc_alloc);
1143 srenew(li->bla,2*li->nc_alloc);
1144 srenew(li->blc,li->nc_alloc);
1145 srenew(li->blc1,li->nc_alloc);
1146 srenew(li->blnr,li->nc_alloc+1);
1147 srenew(li->bllen,li->nc_alloc);
1148 srenew(li->tmpv,li->nc_alloc);
1149 srenew(li->tmp1,li->nc_alloc);
1150 srenew(li->tmp2,li->nc_alloc);
1151 srenew(li->tmp3,li->nc_alloc);
1152 srenew(li->tmp4,li->nc_alloc);
1153 srenew(li->mlambda,li->nc_alloc);
1154 if (li->ncg_triangle > 0)
1156 /* This is allocating too much, but it is difficult to improve */
1157 srenew(li->triangle,li->nc_alloc);
1158 srenew(li->tri_bits,li->nc_alloc);
1162 iatom = idef->il[F_CONSTR].iatoms;
1164 ncc_alloc = li->ncc_alloc;
1167 ni = idef->il[F_CONSTR].nr/3;
1171 li->blnr[con] = nconnect;
1178 lenA = idef->iparams[type].constr.dA;
1179 lenB = idef->iparams[type].constr.dB;
1180 /* Skip the flexible constraints when not doing dynamics */
1181 if (bDynamics || lenA!=0 || lenB!=0)
1183 li->bllen0[con] = lenA;
1184 li->ddist[con] = lenB - lenA;
1185 /* Set the length to the topology A length */
1186 li->bllen[con] = li->bllen0[con];
1187 li->bla[2*con] = a1;
1188 li->bla[2*con+1] = a2;
1189 /* Construct the constraint connection matrix blbnb */
1190 for(k=at2con.index[a1-start]; k<at2con.index[a1-start+1]; k++)
1192 concon = at2con.a[k];
1195 if (nconnect >= ncc_alloc)
1197 ncc_alloc = over_alloc_small(nconnect+1);
1198 srenew(li->blbnb,ncc_alloc);
1200 li->blbnb[nconnect++] = concon;
1203 for(k=at2con.index[a2-start]; k<at2con.index[a2-start+1]; k++)
1205 concon = at2con.a[k];
1208 if (nconnect+1 > ncc_alloc)
1210 ncc_alloc = over_alloc_small(nconnect+1);
1211 srenew(li->blbnb,ncc_alloc);
1213 li->blbnb[nconnect++] = concon;
1216 li->blnr[con+1] = nconnect;
1220 /* Order the blbnb matrix to optimize memory access */
1221 qsort(&(li->blbnb[li->blnr[con]]),li->blnr[con+1]-li->blnr[con],
1222 sizeof(li->blbnb[0]),int_comp);
1224 /* Increase the constraint count */
1229 done_blocka(&at2con);
1231 /* This is the real number of constraints,
1232 * without dynamics the flexible constraints are not present.
1236 li->ncc = li->blnr[con];
1239 /* Since the matrix is static, we can free some memory */
1240 ncc_alloc = li->ncc;
1241 srenew(li->blbnb,ncc_alloc);
1244 if (ncc_alloc > li->ncc_alloc)
1246 li->ncc_alloc = ncc_alloc;
1247 srenew(li->blmf,li->ncc_alloc);
1248 srenew(li->blmf1,li->ncc_alloc);
1249 srenew(li->tmpncc,li->ncc_alloc);
1254 fprintf(debug,"Number of constraints is %d, couplings %d\n",
1261 li->th[0].b1 = li->nc;
1265 lincs_thread_setup(li,md->nr);
1268 set_lincs_matrix(li,md->invmass,md->lambda);
1271 static void lincs_warning(FILE *fplog,
1272 gmx_domdec_t *dd,rvec *x,rvec *xprime,t_pbc *pbc,
1273 int ncons,int *bla,real *bllen,real wangle,
1274 int maxwarn,int *warncount)
1278 real wfac,d0,d1,cosine;
1281 wfac=cos(DEG2RAD*wangle);
1283 sprintf(buf,"bonds that rotated more than %g degrees:\n"
1284 " atom 1 atom 2 angle previous, current, constraint length\n",
1286 fprintf(stderr,"%s",buf);
1289 fprintf(fplog,"%s",buf);
1292 for(b=0;b<ncons;b++)
1298 pbc_dx_aiuc(pbc,x[i],x[j],v0);
1299 pbc_dx_aiuc(pbc,xprime[i],xprime[j],v1);
1303 rvec_sub(x[i],x[j],v0);
1304 rvec_sub(xprime[i],xprime[j],v1);
1308 cosine = iprod(v0,v1)/(d0*d1);
1311 sprintf(buf," %6d %6d %5.1f %8.4f %8.4f %8.4f\n",
1312 ddglatnr(dd,i),ddglatnr(dd,j),
1313 RAD2DEG*acos(cosine),d0,d1,bllen[b]);
1314 fprintf(stderr,"%s",buf);
1317 fprintf(fplog,"%s",buf);
1319 if (!gmx_isfinite(d1))
1321 gmx_fatal(FARGS,"Bond length not finite.");
1327 if (*warncount > maxwarn)
1329 too_many_constraint_warnings(econtLINCS,*warncount);
1333 static void cconerr(gmx_domdec_t *dd,
1334 int ncons,int *bla,real *bllen,rvec *x,t_pbc *pbc,
1335 real *ncons_loc,real *ssd,real *max,int *imax)
1337 real len,d,ma,ssd2,r2;
1338 int *nlocat,count,b,im;
1341 if (dd && dd->constraints)
1343 nlocat = dd_constraints_nlocalatoms(dd);
1354 for(b=0;b<ncons;b++)
1358 pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);
1361 rvec_sub(x[bla[2*b]],x[bla[2*b+1]],dx);
1364 len = r2*gmx_invsqrt(r2);
1365 d = fabs(len/bllen[b]-1);
1366 if (d > ma && (nlocat==NULL || nlocat[b]))
1378 ssd2 += nlocat[b]*d*d;
1383 *ncons_loc = (nlocat ? 0.5 : 1)*count;
1384 *ssd = (nlocat ? 0.5 : 1)*ssd2;
1389 static void dump_conf(gmx_domdec_t *dd,struct gmx_lincsdata *li,
1391 char *name,gmx_bool bAll,rvec *x,matrix box)
1397 dd_get_constraint_range(dd,&ac0,&ac1);
1399 sprintf(str,"%s_%d_%d_%d.pdb",name,dd->ci[XX],dd->ci[YY],dd->ci[ZZ]);
1400 fp = gmx_fio_fopen(str,"w");
1401 fprintf(fp,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1\n",
1402 10*norm(box[XX]),10*norm(box[YY]),10*norm(box[ZZ]),
1404 for(i=0; i<ac1; i++)
1406 if (i < dd->nat_home || (bAll && i >= ac0 && i < ac1))
1408 fprintf(fp,"%-6s%5u %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
1409 "ATOM",ddglatnr(dd,i),"C","ALA",' ',i+1,
1410 10*x[i][XX],10*x[i][YY],10*x[i][ZZ],
1411 1.0,i<dd->nat_tot ? 0.0 : 1.0);
1416 for(i=0; i<li->nc; i++)
1418 fprintf(fp,"CONECT%5d%5d\n",
1419 ddglatnr(dd,li->bla[2*i]),
1420 ddglatnr(dd,li->bla[2*i+1]));
1426 gmx_bool constrain_lincs(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
1428 gmx_large_int_t step,
1429 struct gmx_lincsdata *lincsd,t_mdatoms *md,
1431 rvec *x,rvec *xprime,rvec *min_proj,
1432 matrix box,t_pbc *pbc,
1433 real lambda,real *dvdlambda,
1435 gmx_bool bCalcVir,tensor vir_r_m_dr,
1438 int maxwarn,int *warncount)
1440 char buf[STRLEN],buf2[22],buf3[STRLEN];
1441 int i,warn,p_imax,error;
1442 real ncons_loc,p_ssd,p_max=0;
1448 if (lincsd->nc == 0 && cr->dd == NULL)
1452 lincsd->rmsd_data[0] = 0;
1453 if (ir->eI == eiSD2 && v == NULL)
1461 lincsd->rmsd_data[i] = 0;
1467 if (econq == econqCoord)
1469 if (ir->efep != efepNO)
1471 if (md->nMassPerturbed && lincsd->matlam != md->lambda)
1473 set_lincs_matrix(lincsd,md->invmass,md->lambda);
1476 for(i=0; i<lincsd->nc; i++)
1478 lincsd->bllen[i] = lincsd->bllen0[i] + lambda*lincsd->ddist[i];
1482 if (lincsd->ncg_flex)
1484 /* Set the flexible constraint lengths to the old lengths */
1487 for(i=0; i<lincsd->nc; i++)
1489 if (lincsd->bllen[i] == 0) {
1490 pbc_dx_aiuc(pbc,x[lincsd->bla[2*i]],x[lincsd->bla[2*i+1]],dx);
1491 lincsd->bllen[i] = norm(dx);
1497 for(i=0; i<lincsd->nc; i++)
1499 if (lincsd->bllen[i] == 0)
1502 sqrt(distance2(x[lincsd->bla[2*i]],
1503 x[lincsd->bla[2*i+1]]));
1511 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc,
1512 &ncons_loc,&p_ssd,&p_max,&p_imax);
1515 /* This warn var can be updated by multiple threads
1516 * at the same time. But as we only need to detect
1517 * if a warning occured or not, this is not an issue.
1521 /* The OpenMP parallel region of constrain_lincs for coords */
1522 #pragma omp parallel num_threads(lincsd->nth)
1524 int th=gmx_omp_get_thread_num();
1526 clear_mat(lincsd->th[th].vir_r_m_dr);
1528 do_lincs(x,xprime,box,pbc,lincsd,th,
1530 bCalcVir || (ir->efep != efepNO),
1531 ir->LincsWarnAngle,&warn,
1533 th==0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
1536 if (ir->efep != efepNO)
1540 dt_2 = 1.0/(ir->delta_t*ir->delta_t);
1541 for(i=0; (i<lincsd->nc); i++)
1543 dvdl -= lincsd->mlambda[i]*dt_2*lincsd->ddist[i];
1548 if (bLog && fplog && lincsd->nc > 0)
1550 fprintf(fplog," Rel. Constraint Deviation: RMS MAX between atoms\n");
1551 fprintf(fplog," Before LINCS %.6f %.6f %6d %6d\n",
1552 sqrt(p_ssd/ncons_loc),p_max,
1553 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1554 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1558 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc,
1559 &ncons_loc,&p_ssd,&p_max,&p_imax);
1560 /* Check if we are doing the second part of SD */
1561 if (ir->eI == eiSD2 && v == NULL)
1569 lincsd->rmsd_data[0] = ncons_loc;
1570 lincsd->rmsd_data[i] = p_ssd;
1574 lincsd->rmsd_data[0] = 0;
1575 lincsd->rmsd_data[1] = 0;
1576 lincsd->rmsd_data[2] = 0;
1578 if (bLog && fplog && lincsd->nc > 0)
1581 " After LINCS %.6f %.6f %6d %6d\n\n",
1582 sqrt(p_ssd/ncons_loc),p_max,
1583 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1584 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1591 cconerr(cr->dd,lincsd->nc,lincsd->bla,lincsd->bllen,xprime,pbc,
1592 &ncons_loc,&p_ssd,&p_max,&p_imax);
1595 sprintf(buf3," in simulation %d", cr->ms->sim);
1601 sprintf(buf,"\nStep %s, time %g (ps) LINCS WARNING%s\n"
1602 "relative constraint deviation after LINCS:\n"
1603 "rms %.6f, max %.6f (between atoms %d and %d)\n",
1604 gmx_step_str(step,buf2),ir->init_t+step*ir->delta_t,
1606 sqrt(p_ssd/ncons_loc),p_max,
1607 ddglatnr(cr->dd,lincsd->bla[2*p_imax]),
1608 ddglatnr(cr->dd,lincsd->bla[2*p_imax+1]));
1611 fprintf(fplog,"%s",buf);
1613 fprintf(stderr,"%s",buf);
1614 lincs_warning(fplog,cr->dd,x,xprime,pbc,
1615 lincsd->nc,lincsd->bla,lincsd->bllen,
1616 ir->LincsWarnAngle,maxwarn,warncount);
1618 bOK = (p_max < 0.5);
1621 if (lincsd->ncg_flex) {
1622 for(i=0; (i<lincsd->nc); i++)
1623 if (lincsd->bllen0[i] == 0 && lincsd->ddist[i] == 0)
1624 lincsd->bllen[i] = 0;
1629 /* The OpenMP parallel region of constrain_lincs for derivatives */
1630 #pragma omp parallel num_threads(lincsd->nth)
1632 int th=gmx_omp_get_thread_num();
1634 do_lincsp(x,xprime,min_proj,pbc,lincsd,th,
1635 md->invmass,econq,ir->efep != efepNO ? dvdlambda : NULL,
1636 bCalcVir,th==0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
1640 if (bCalcVir && lincsd->nth > 1)
1642 for(i=1; i<lincsd->nth; i++)
1644 m_add(vir_r_m_dr,lincsd->th[i].vir_r_m_dr,vir_r_m_dr);
1648 /* count assuming nit=1 */
1649 inc_nrnb(nrnb,eNR_LINCS,lincsd->nc);
1650 inc_nrnb(nrnb,eNR_LINCSMAT,(2+lincsd->nOrder)*lincsd->ncc);
1651 if (lincsd->ntriangle > 0)
1653 inc_nrnb(nrnb,eNR_LINCSMAT,lincsd->nOrder*lincsd->ncc_triangle);
1657 inc_nrnb(nrnb,eNR_CONSTR_V,lincsd->nc*2);
1661 inc_nrnb(nrnb,eNR_CONSTR_VIR,lincsd->nc);