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"
58 int b0
; /* first constraint for this thread */
59 int b1
; /* b1-1 is the last constraint for this thread */
60 int nind
; /* number of indices */
61 int *ind
; /* constraint index for updating atom data */
62 int nind_r
; /* number of indices */
63 int *ind_r
; /* constraint index for updating atom data */
64 int ind_nalloc
; /* allocation size of ind and ind_r */
67 typedef struct gmx_lincsdata
{
68 int ncg
; /* the global number of constraints */
69 int ncg_flex
; /* the global number of flexible constraints */
70 int ncg_triangle
;/* the global number of constraints in triangles */
71 int nIter
; /* the number of iterations */
72 int nOrder
; /* the order of the matrix expansion */
73 int nc
; /* the number of constraints */
74 int nc_alloc
; /* the number we allocated memory for */
75 int ncc
; /* the number of constraint connections */
76 int ncc_alloc
; /* the number we allocated memory for */
77 real matlam
; /* the FE lambda value used for filling blc and blmf */
78 real
*bllen0
; /* the reference distance in topology A */
79 real
*ddist
; /* the reference distance in top B - the r.d. in top A */
80 int *bla
; /* the atom pairs involved in the constraints */
81 real
*blc
; /* 1/sqrt(invmass1 + invmass2) */
82 real
*blc1
; /* as blc, but with all masses 1 */
83 int *blnr
; /* index into blbnb and blmf */
84 int *blbnb
; /* list of constraint connections */
85 int ntriangle
; /* the local number of constraints in triangles */
86 int *triangle
; /* the list of triangle constraints */
87 int *tri_bits
; /* the bits tell if the matrix element should be used */
88 int ncc_triangle
;/* the number of constraint connections in triangles */
89 real
*blmf
; /* matrix of mass factors for constraint connections */
90 real
*blmf1
; /* as blmf, but with all masses 1 */
91 real
*bllen
; /* the reference bond length */
92 int nth
; /* The number of threads doing LINCS */
93 lincs_thread_t
*th
; /* LINCS thread division */
94 unsigned *atf
; /* atom flags for thread parallelization */
95 int atf_nalloc
; /* allocation size of atf */
96 /* arrays for temporary storage in the LINCS algorithm */
103 real
*mlambda
; /* the Lagrange multipliers * -1 */
104 /* storage for the constraint RMS relative deviation output */
108 real
*lincs_rmsd_data(struct gmx_lincsdata
*lincsd
)
110 return lincsd
->rmsd_data
;
113 real
lincs_rmsd(struct gmx_lincsdata
*lincsd
,gmx_bool bSD2
)
115 if (lincsd
->rmsd_data
[0] > 0)
117 return sqrt(lincsd
->rmsd_data
[bSD2
? 2 : 1]/lincsd
->rmsd_data
[0]);
125 /* Do a set of nrec LINCS matrix multiplications.
126 * This function will return with up to date thread-local
127 * constraint data, without an OpenMP barrier.
129 static void lincs_matrix_expand(const struct gmx_lincsdata
*lincsd
,
132 real
*rhs1
,real
*rhs2
,real
*sol
)
134 int nrec
,rec
,b
,j
,n
,nr0
,nr1
;
136 int ntriangle
,tb
,bits
;
137 const int *blnr
=lincsd
->blnr
,*blbnb
=lincsd
->blbnb
;
138 const int *triangle
=lincsd
->triangle
,*tri_bits
=lincsd
->tri_bits
;
140 ntriangle
= lincsd
->ntriangle
;
141 nrec
= lincsd
->nOrder
;
143 for(rec
=0; rec
<nrec
; rec
++)
149 for(n
=blnr
[b
]; n
<blnr
[b
+1]; n
++)
152 mvb
= mvb
+ blcc
[n
]*rhs1
[j
];
155 sol
[b
] = sol
[b
] + mvb
;
160 } /* nrec*(ncons+2*nrtot) flops */
164 /* Perform an extra nrec recursions for only the constraints
165 * involved in rigid triangles.
166 * In this way their accuracy should come close to those of the other
167 * constraints, since traingles of constraints can produce eigenvalues
168 * around 0.7, while the effective eigenvalue for bond constraints
169 * is around 0.4 (and 0.7*0.7=0.5).
171 /* We need to copy the temporary array, since only the elements
172 * for constraints involved in triangles are updated and then
173 * the pointers are swapped. This saving copying the whole arrary.
174 * We need barrier as other threads might still be reading from rhs2.
184 for(rec
=0; rec
<nrec
; rec
++)
186 for(tb
=0; tb
<ntriangle
; tb
++)
193 for(n
=nr0
; n
<nr1
; n
++)
195 if (bits
& (1<<(n
-nr0
)))
198 mvb
= mvb
+ blcc
[n
]*rhs1
[j
];
202 sol
[b
] = sol
[b
] + mvb
;
208 } /* flops count is missing here */
210 /* We need a barrier here as the calling routine will continue
211 * to operate on the thread-local constraints without barrier.
217 static void lincs_update_atoms_noind(int ncons
,const int *bla
,
219 const real
*fac
,rvec
*r
,
224 real mvb
,im1
,im2
,tmp0
,tmp1
,tmp2
;
226 for(b
=0; b
<ncons
; b
++)
242 } /* 16 ncons flops */
245 static void lincs_update_atoms_ind(int ncons
,const int *ind
,const int *bla
,
247 const real
*fac
,rvec
*r
,
252 real mvb
,im1
,im2
,tmp0
,tmp1
,tmp2
;
254 for(bi
=0; bi
<ncons
; bi
++)
271 } /* 16 ncons flops */
274 static void lincs_update_atoms(struct gmx_lincsdata
*li
,int th
,
276 const real
*fac
,rvec
*r
,
282 /* Single thread, we simply update for all constraints */
283 lincs_update_atoms_noind(li
->nc
,li
->bla
,prefac
,fac
,r
,invmass
,x
);
287 /* Update the atom vector components for our thread local
288 * constraints that only access our local atom range.
289 * This can be done without a barrier.
291 lincs_update_atoms_ind(li
->th
[th
].nind
,li
->th
[th
].ind
,
292 li
->bla
,prefac
,fac
,r
,invmass
,x
);
294 if (li
->th
[li
->nth
].nind
> 0)
296 /* Update the constraints that operate on atoms
297 * in multiple thread atom blocks on the master thread.
302 lincs_update_atoms_ind(li
->th
[li
->nth
].nind
,
304 li
->bla
,prefac
,fac
,r
,invmass
,x
);
310 static void do_lincsp(rvec
*x
,rvec
*f
,rvec
*fp
,t_pbc
*pbc
,
311 struct gmx_lincsdata
*lincsd
,real
*invmass
,
312 int econq
,real
*dvdlambda
,
313 gmx_bool bCalcVir
,tensor rmdf
)
316 real tmp0
,tmp1
,tmp2
,im1
,im2
,mvb
,rlen
,len
,wfac
,lam
;
318 int ncons
,*bla
,*blnr
,*blbnb
;
320 real
*blc
,*blmf
,*blcc
,*rhs1
,*rhs2
,*sol
;
326 blbnb
= lincsd
->blbnb
;
327 if (econq
!= econqForce
)
329 /* Use mass-weighted parameters */
335 /* Use non mass-weighted parameters */
337 blmf
= lincsd
->blmf1
;
339 blcc
= lincsd
->tmpncc
;
344 if (econq
!= econqForce
)
349 /* Compute normalized i-j vectors */
352 for(b
=0; b
<ncons
; b
++)
354 pbc_dx_aiuc(pbc
,x
[bla
[2*b
]],x
[bla
[2*b
+1]],dx
);
360 for(b
=0; b
<ncons
; b
++)
362 rvec_sub(x
[bla
[2*b
]],x
[bla
[2*b
+1]],dx
);
364 } /* 16 ncons flops */
367 for(b
=0; b
<ncons
; b
++)
374 for(n
=blnr
[b
]; n
<blnr
[b
+1]; n
++)
377 blcc
[n
] = blmf
[n
]*(tmp0
*r
[k
][0] + tmp1
*r
[k
][1] + tmp2
*r
[k
][2]);
379 mvb
= blc
[b
]*(tmp0
*(f
[i
][0] - f
[j
][0]) +
380 tmp1
*(f
[i
][1] - f
[j
][1]) +
381 tmp2
*(f
[i
][2] - f
[j
][2]));
386 /* Together: 23*ncons + 6*nrtot flops */
388 lincs_matrix_expand(lincsd
,0,ncons
,blcc
,rhs1
,rhs2
,sol
);
389 /* nrec*(ncons+2*nrtot) flops */
391 if (econq
!= econqForce
)
393 for(b
=0; b
<ncons
; b
++)
395 /* With econqDeriv_FlexCon only use the flexible constraints */
396 if (econq
!= econqDeriv_FlexCon
||
397 (lincsd
->bllen0
[b
] == 0 && lincsd
->ddist
[b
] == 0))
407 fp
[i
][0] -= tmp0
*im1
;
408 fp
[i
][1] -= tmp1
*im1
;
409 fp
[i
][2] -= tmp2
*im1
;
410 fp
[j
][0] += tmp0
*im2
;
411 fp
[j
][1] += tmp1
*im2
;
412 fp
[j
][2] += tmp2
*im2
;
415 /* This is only correct with forces and invmass=1 */
416 *dvdlambda
-= mvb
*lincsd
->ddist
[b
];
419 } /* 16 ncons flops */
423 for(b
=0; b
<ncons
; b
++)
439 *dvdlambda
-= mvb
*lincsd
->ddist
[b
];
447 /* Constraint virial,
448 * determines sum r_bond x delta f,
449 * where delta f is the constraint correction
450 * of the quantity that is being constrained.
452 for(b
=0; b
<ncons
; b
++)
454 mvb
= lincsd
->bllen
[b
]*blc
[b
]*sol
[b
];
460 rmdf
[i
][j
] += tmp1
*r
[b
][j
];
463 } /* 23 ncons flops */
467 static void do_lincs(rvec
*x
,rvec
*xp
,matrix box
,t_pbc
*pbc
,
468 struct gmx_lincsdata
*lincsd
,int th
,
471 gmx_bool bCalcLambda
,
472 real wangle
,int *warn
,
474 gmx_bool bCalcVir
,tensor rmdr
)
476 int b0
,b1
,b
,i
,j
,k
,n
,iter
;
477 real tmp0
,tmp1
,tmp2
,im1
,im2
,mvb
,rlen
,len
,len2
,dlen2
,wfac
;
479 int ncons
,*bla
,*blnr
,*blbnb
;
481 real
*blc
,*blmf
,*bllen
,*blcc
,*rhs1
,*rhs2
,*sol
,*blc_sol
,*mlambda
;
484 b0
= lincsd
->th
[th
].b0
;
485 b1
= lincsd
->th
[th
].b1
;
491 blbnb
= lincsd
->blbnb
;
494 bllen
= lincsd
->bllen
;
495 blcc
= lincsd
->tmpncc
;
499 blc_sol
= lincsd
->tmp4
;
500 mlambda
= lincsd
->mlambda
;
502 if (DOMAINDECOMP(cr
) && cr
->dd
->constraints
)
504 nlocat
= dd_constraints_nlocalatoms(cr
->dd
);
506 else if (PARTDECOMP(cr
))
508 nlocat
= pd_constraints_nlocalatoms(cr
->pd
);
519 /* Compute normalized i-j vectors */
522 pbc_dx_aiuc(pbc
,x
[bla
[2*b
]],x
[bla
[2*b
+1]],dx
);
528 for(n
=blnr
[b
]; n
<blnr
[b
+1]; n
++)
530 blcc
[n
] = blmf
[n
]*iprod(r
[b
],r
[blbnb
[n
]]);
532 pbc_dx_aiuc(pbc
,xp
[bla
[2*b
]],xp
[bla
[2*b
+1]],dx
);
533 mvb
= blc
[b
]*(iprod(r
[b
],dx
) - bllen
[b
]);
540 /* Compute normalized i-j vectors */
545 tmp0
= x
[i
][0] - x
[j
][0];
546 tmp1
= x
[i
][1] - x
[j
][1];
547 tmp2
= x
[i
][2] - x
[j
][2];
548 rlen
= gmx_invsqrt(tmp0
*tmp0
+tmp1
*tmp1
+tmp2
*tmp2
);
552 } /* 16 ncons flops */
563 for(n
=blnr
[b
]; n
<blnr
[b
+1]; n
++)
566 blcc
[n
] = blmf
[n
]*(tmp0
*r
[k
][0] + tmp1
*r
[k
][1] + tmp2
*r
[k
][2]);
568 mvb
= blc
[b
]*(tmp0
*(xp
[i
][0] - xp
[j
][0]) +
569 tmp1
*(xp
[i
][1] - xp
[j
][1]) +
570 tmp2
*(xp
[i
][2] - xp
[j
][2]) - len
);
575 /* Together: 26*ncons + 6*nrtot flops */
578 lincs_matrix_expand(lincsd
,b0
,b1
,blcc
,rhs1
,rhs2
,sol
);
579 /* nrec*(ncons+2*nrtot) flops */
583 mlambda
[b
] = blc
[b
]*sol
[b
];
586 /* Update the coordinates */
587 lincs_update_atoms(lincsd
,th
,1.0,mlambda
,r
,invmass
,xp
);
590 ******** Correction for centripetal effects ********
593 wfac
= cos(DEG2RAD
*wangle
);
596 for(iter
=0; iter
<lincsd
->nIter
; iter
++)
598 if ((DOMAINDECOMP(cr
) && cr
->dd
->constraints
) ||
604 /* Communicate the corrected non-local coordinates */
605 if (DOMAINDECOMP(cr
))
607 dd_move_x_constraints(cr
->dd
,box
,xp
,NULL
);
611 pd_move_x_constraints(cr
,xp
,NULL
);
622 pbc_dx_aiuc(pbc
,xp
[bla
[2*b
]],xp
[bla
[2*b
+1]],dx
);
626 rvec_sub(xp
[bla
[2*b
]],xp
[bla
[2*b
+1]],dx
);
629 dlen2
= 2*len2
- norm2(dx
);
630 if (dlen2
< wfac
*len2
&& (nlocat
==NULL
|| nlocat
[b
]))
636 mvb
= blc
[b
]*(len
- dlen2
*gmx_invsqrt(dlen2
));
644 } /* 20*ncons flops */
646 lincs_matrix_expand(lincsd
,b0
,b1
,blcc
,rhs1
,rhs2
,sol
);
647 /* nrec*(ncons+2*nrtot) flops */
656 /* Update the coordinates */
657 lincs_update_atoms(lincsd
,th
,1.0,blc_sol
,r
,invmass
,xp
);
659 /* nit*ncons*(37+9*nrec) flops */
663 /* Update the velocities */
664 lincs_update_atoms(lincsd
,th
,invdt
,mlambda
,r
,invmass
,v
);
668 if (nlocat
&& bCalcLambda
)
670 /* Only account for local atoms */
673 mlambda
[b
] *= 0.5*nlocat
[b
];
682 /* Constraint virial */
683 for(b
=0; b
<ncons
; b
++)
685 tmp0
= -bllen
[b
]*mlambda
[b
];
691 rmdr
[i
][j
] -= tmp1
*r
[b
][j
];
694 } /* 22 ncons flops */
699 * 26*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
700 * + nit * (20*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
702 * (26+nrec)*ncons + (6+2*nrec)*nrtot
703 * + nit * ((37+nrec)*ncons + 2*nrec*nrtot)
705 * (63+nrec)*ncons + (6+4*nrec)*nrtot
709 void set_lincs_matrix(struct gmx_lincsdata
*li
,real
*invmass
,real lambda
)
711 int i
,a1
,a2
,n
,k
,sign
,center
;
713 const real invsqrt2
=0.7071067811865475244;
715 for(i
=0; (i
<li
->nc
); i
++)
719 li
->blc
[i
] = gmx_invsqrt(invmass
[a1
] + invmass
[a2
]);
720 li
->blc1
[i
] = invsqrt2
;
723 /* Construct the coupling coefficient matrix blmf */
725 li
->ncc_triangle
= 0;
726 for(i
=0; (i
<li
->nc
); i
++)
730 for(n
=li
->blnr
[i
]; (n
<li
->blnr
[i
+1]); n
++)
733 if (a1
== li
->bla
[2*k
] || a2
== li
->bla
[2*k
+1])
741 if (a1
== li
->bla
[2*k
] || a1
== li
->bla
[2*k
+1])
751 li
->blmf
[n
] = sign
*invmass
[center
]*li
->blc
[i
]*li
->blc
[k
];
752 li
->blmf1
[n
] = sign
*0.5;
753 if (li
->ncg_triangle
> 0)
755 /* Look for constraint triangles */
756 for(nk
=li
->blnr
[k
]; (nk
<li
->blnr
[k
+1]); nk
++)
759 if (kk
!= i
&& kk
!= k
&&
760 (li
->bla
[2*kk
] == end
|| li
->bla
[2*kk
+1] == end
))
762 if (li
->ntriangle
== 0 ||
763 li
->triangle
[li
->ntriangle
-1] < i
)
765 /* Add this constraint to the triangle list */
766 li
->triangle
[li
->ntriangle
] = i
;
767 li
->tri_bits
[li
->ntriangle
] = 0;
769 if (li
->blnr
[i
+1] - li
->blnr
[i
] > sizeof(li
->tri_bits
[0])*8 - 1)
771 gmx_fatal(FARGS
,"A constraint is connected to %d constraints, this is more than the %d allowed for constraints participating in triangles",
772 li
->blnr
[i
+1] - li
->blnr
[i
],
773 sizeof(li
->tri_bits
[0])*8-1);
776 li
->tri_bits
[li
->ntriangle
-1] |= (1<<(n
-li
->blnr
[i
]));
786 fprintf(debug
,"Of the %d constraints %d participate in triangles\n",
787 li
->nc
,li
->ntriangle
);
788 fprintf(debug
,"There are %d couplings of which %d in triangles\n",
789 li
->ncc
,li
->ncc_triangle
);
793 * so we know with which lambda value the masses have been set.
798 static int count_triangle_constraints(t_ilist
*ilist
,t_blocka
*at2con
)
801 int c0
,a00
,a01
,n1
,c1
,a10
,a11
,ac1
,n2
,c2
,a20
,a21
;
804 t_iatom
*ia1
,*ia2
,*iap
;
806 ncon1
= ilist
[F_CONSTR
].nr
/3;
807 ncon_tot
= ncon1
+ ilist
[F_CONSTRNC
].nr
/3;
809 ia1
= ilist
[F_CONSTR
].iatoms
;
810 ia2
= ilist
[F_CONSTRNC
].iatoms
;
813 for(c0
=0; c0
<ncon_tot
; c0
++)
816 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,c0
);
819 for(n1
=at2con
->index
[a01
]; n1
<at2con
->index
[a01
+1]; n1
++)
824 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,c1
);
835 for(n2
=at2con
->index
[ac1
]; n2
<at2con
->index
[ac1
+1]; n2
++)
838 if (c2
!= c0
&& c2
!= c1
)
840 iap
= constr_iatomptr(ncon1
,ia1
,ia2
,c2
);
843 if (a20
== a00
|| a21
== a00
)
857 return ncon_triangle
;
860 static int int_comp(const void *a
,const void *b
)
862 return (*(int *)a
) - (*(int *)b
);
865 gmx_lincsdata_t
init_lincs(FILE *fplog
,gmx_mtop_t
*mtop
,
866 int nflexcon_global
,t_blocka
*at2con
,
867 gmx_bool bPLINCS
,int nIter
,int nProjOrder
)
869 struct gmx_lincsdata
*li
;
875 fprintf(fplog
,"\nInitializing%s LINear Constraint Solver\n",
876 bPLINCS
? " Parallel" : "");
882 gmx_mtop_ftype_count(mtop
,F_CONSTR
) +
883 gmx_mtop_ftype_count(mtop
,F_CONSTRNC
);
884 li
->ncg_flex
= nflexcon_global
;
886 li
->ncg_triangle
= 0;
887 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
889 molt
= &mtop
->moltype
[mtop
->molblock
[mb
].type
];
891 mtop
->molblock
[mb
].nmol
*
892 count_triangle_constraints(molt
->ilist
,
893 &at2con
[mtop
->molblock
[mb
].type
]);
897 li
->nOrder
= nProjOrder
;
899 /* LINCS can run on any number of threads.
900 * Currently the number is fixed for the whole simulation,
901 * but it could be set in set_lincs().
903 li
->nth
= gmx_omp_nthreads_get(emntLINCS
);
910 /* Allocate an extra elements for "thread-overlap" constraints */
911 snew(li
->th
,li
->nth
+1);
915 fprintf(debug
,"LINCS: using %d threads\n",li
->nth
);
918 if (bPLINCS
|| li
->ncg_triangle
> 0)
920 please_cite(fplog
,"Hess2008a");
924 please_cite(fplog
,"Hess97a");
929 fprintf(fplog
,"The number of constraints is %d\n",li
->ncg
);
932 fprintf(fplog
,"There are inter charge-group constraints,\n"
933 "will communicate selected coordinates each lincs iteration\n");
935 if (li
->ncg_triangle
> 0)
938 "%d constraints are involved in constraint triangles,\n"
939 "will apply an additional matrix expansion of order %d for couplings\n"
940 "between constraints inside triangles\n",
941 li
->ncg_triangle
,li
->nOrder
);
948 /* Sets up the work division over the threads */
949 static void lincs_thread_setup(struct gmx_lincsdata
*li
,int natoms
)
951 lincs_thread_t
*li_m
;
956 if (natoms
> li
->atf_nalloc
)
958 li
->atf_nalloc
= over_alloc_large(natoms
);
959 srenew(li
->atf
,li
->atf_nalloc
);
963 /* Clear the atom flags */
964 for(a
=0; a
<natoms
; a
++)
969 for(th
=0; th
<li
->nth
; th
++)
971 lincs_thread_t
*li_th
;
976 /* The constraints are divided equally over the threads */
977 li_th
->b0
= (li
->nc
* th
)/li
->nth
;
978 li_th
->b1
= (li
->nc
*(th
+1))/li
->nth
;
980 if (th
< sizeof(*atf
)*8)
982 /* For each atom set a flag for constraints from each */
983 for(b
=li_th
->b0
; b
<li_th
->b1
; b
++)
985 atf
[li
->bla
[b
*2] ] |= (1U<<th
);
986 atf
[li
->bla
[b
*2+1]] |= (1U<<th
);
991 #pragma omp parallel for num_threads(li->nth) schedule(static)
992 for(th
=0; th
<li
->nth
; th
++)
994 lincs_thread_t
*li_th
;
1000 if (li_th
->b1
- li_th
->b0
> li_th
->ind_nalloc
)
1002 li_th
->ind_nalloc
= over_alloc_large(li_th
->b1
-li_th
->b0
);
1003 srenew(li_th
->ind
,li_th
->ind_nalloc
);
1004 srenew(li_th
->ind_r
,li_th
->ind_nalloc
);
1007 if (th
< sizeof(*atf
)*8)
1009 mask
= (1U<<th
) - 1U;
1013 for(b
=li_th
->b0
; b
<li_th
->b1
; b
++)
1015 /* We let the constraint with the lowest thread index
1016 * operate on atoms with constraints from multiple threads.
1018 if (((atf
[li
->bla
[b
*2]] & mask
) == 0) &&
1019 ((atf
[li
->bla
[b
*2+1]] & mask
) == 0))
1021 /* Add the constraint to the local atom update index */
1022 li_th
->ind
[li_th
->nind
++] = b
;
1026 /* Add the constraint to the rest block */
1027 li_th
->ind_r
[li_th
->nind_r
++] = b
;
1033 /* We are out of bits, assign all constraints to rest */
1034 for(b
=li_th
->b0
; b
<li_th
->b1
; b
++)
1036 li_th
->ind_r
[li_th
->nind_r
++] = b
;
1041 /* We need to copy all constraints which have not be assigned
1042 * to a thread to a separate list which will be handled by one thread.
1044 li_m
= &li
->th
[li
->nth
];
1047 for(th
=0; th
<li
->nth
; th
++)
1049 lincs_thread_t
*li_th
;
1052 li_th
= &li
->th
[th
];
1054 if (li_m
->nind
+ li_th
->nind_r
> li_m
->ind_nalloc
)
1056 li_m
->ind_nalloc
= over_alloc_large(li_m
->nind
+li_th
->nind_r
);
1057 srenew(li_m
->ind
,li_m
->ind_nalloc
);
1060 for(b
=0; b
<li_th
->nind_r
; b
++)
1062 li_m
->ind
[li_m
->nind
++] = li_th
->ind_r
[b
];
1067 fprintf(debug
,"LINCS thread %d: %d constraints\n",
1074 fprintf(debug
,"LINCS thread r: %d constraints\n",
1080 void set_lincs(t_idef
*idef
,t_mdatoms
*md
,
1081 gmx_bool bDynamics
,t_commrec
*cr
,
1082 struct gmx_lincsdata
*li
)
1084 int start
,natoms
,nflexcon
;
1087 int i
,k
,ncc_alloc
,ni
,con
,nconnect
,concon
;
1094 /* Zero the thread index ranges.
1095 * Otherwise without local constraints we could return with old ranges.
1097 for(i
=0; i
<li
->nth
; i
++)
1105 li
->th
[li
->nth
].nind
= 0;
1108 /* This is the local topology, so there are only F_CONSTR constraints */
1109 if (idef
->il
[F_CONSTR
].nr
== 0)
1111 /* There are no constraints,
1112 * we do not need to fill any data structures.
1119 fprintf(debug
,"Building the LINCS connectivity\n");
1122 if (DOMAINDECOMP(cr
))
1124 if (cr
->dd
->constraints
)
1126 dd_get_constraint_range(cr
->dd
,&start
,&natoms
);
1130 natoms
= cr
->dd
->nat_home
;
1134 else if(PARTDECOMP(cr
))
1136 pd_get_constraint_range(cr
->pd
,&start
,&natoms
);
1141 natoms
= md
->homenr
;
1143 at2con
= make_at2con(start
,natoms
,idef
->il
,idef
->iparams
,bDynamics
,
1147 if (idef
->il
[F_CONSTR
].nr
/3 > li
->nc_alloc
|| li
->nc_alloc
== 0)
1149 li
->nc_alloc
= over_alloc_dd(idef
->il
[F_CONSTR
].nr
/3);
1150 srenew(li
->bllen0
,li
->nc_alloc
);
1151 srenew(li
->ddist
,li
->nc_alloc
);
1152 srenew(li
->bla
,2*li
->nc_alloc
);
1153 srenew(li
->blc
,li
->nc_alloc
);
1154 srenew(li
->blc1
,li
->nc_alloc
);
1155 srenew(li
->blnr
,li
->nc_alloc
+1);
1156 srenew(li
->bllen
,li
->nc_alloc
);
1157 srenew(li
->tmpv
,li
->nc_alloc
);
1158 srenew(li
->tmp1
,li
->nc_alloc
);
1159 srenew(li
->tmp2
,li
->nc_alloc
);
1160 srenew(li
->tmp3
,li
->nc_alloc
);
1161 srenew(li
->tmp4
,li
->nc_alloc
);
1162 srenew(li
->mlambda
,li
->nc_alloc
);
1163 if (li
->ncg_triangle
> 0)
1165 /* This is allocating too much, but it is difficult to improve */
1166 srenew(li
->triangle
,li
->nc_alloc
);
1167 srenew(li
->tri_bits
,li
->nc_alloc
);
1171 iatom
= idef
->il
[F_CONSTR
].iatoms
;
1173 ncc_alloc
= li
->ncc_alloc
;
1176 ni
= idef
->il
[F_CONSTR
].nr
/3;
1180 li
->blnr
[con
] = nconnect
;
1187 lenA
= idef
->iparams
[type
].constr
.dA
;
1188 lenB
= idef
->iparams
[type
].constr
.dB
;
1189 /* Skip the flexible constraints when not doing dynamics */
1190 if (bDynamics
|| lenA
!=0 || lenB
!=0)
1192 li
->bllen0
[con
] = lenA
;
1193 li
->ddist
[con
] = lenB
- lenA
;
1194 /* Set the length to the topology A length */
1195 li
->bllen
[con
] = li
->bllen0
[con
];
1196 li
->bla
[2*con
] = a1
;
1197 li
->bla
[2*con
+1] = a2
;
1198 /* Construct the constraint connection matrix blbnb */
1199 for(k
=at2con
.index
[a1
-start
]; k
<at2con
.index
[a1
-start
+1]; k
++)
1201 concon
= at2con
.a
[k
];
1204 if (nconnect
>= ncc_alloc
)
1206 ncc_alloc
= over_alloc_small(nconnect
+1);
1207 srenew(li
->blbnb
,ncc_alloc
);
1209 li
->blbnb
[nconnect
++] = concon
;
1212 for(k
=at2con
.index
[a2
-start
]; k
<at2con
.index
[a2
-start
+1]; k
++)
1214 concon
= at2con
.a
[k
];
1217 if (nconnect
+1 > ncc_alloc
)
1219 ncc_alloc
= over_alloc_small(nconnect
+1);
1220 srenew(li
->blbnb
,ncc_alloc
);
1222 li
->blbnb
[nconnect
++] = concon
;
1225 li
->blnr
[con
+1] = nconnect
;
1229 /* Order the blbnb matrix to optimize memory access */
1230 qsort(&(li
->blbnb
[li
->blnr
[con
]]),li
->blnr
[con
+1]-li
->blnr
[con
],
1231 sizeof(li
->blbnb
[0]),int_comp
);
1233 /* Increase the constraint count */
1238 done_blocka(&at2con
);
1240 /* This is the real number of constraints,
1241 * without dynamics the flexible constraints are not present.
1245 li
->ncc
= li
->blnr
[con
];
1248 /* Since the matrix is static, we can free some memory */
1249 ncc_alloc
= li
->ncc
;
1250 srenew(li
->blbnb
,ncc_alloc
);
1253 if (ncc_alloc
> li
->ncc_alloc
)
1255 li
->ncc_alloc
= ncc_alloc
;
1256 srenew(li
->blmf
,li
->ncc_alloc
);
1257 srenew(li
->blmf1
,li
->ncc_alloc
);
1258 srenew(li
->tmpncc
,li
->ncc_alloc
);
1263 fprintf(debug
,"Number of constraints is %d, couplings %d\n",
1270 li
->th
[0].b1
= li
->nc
;
1274 lincs_thread_setup(li
,md
->nr
);
1277 set_lincs_matrix(li
,md
->invmass
,md
->lambda
);
1280 static void lincs_warning(FILE *fplog
,
1281 gmx_domdec_t
*dd
,rvec
*x
,rvec
*xprime
,t_pbc
*pbc
,
1282 int ncons
,int *bla
,real
*bllen
,real wangle
,
1283 int maxwarn
,int *warncount
)
1287 real wfac
,d0
,d1
,cosine
;
1290 wfac
=cos(DEG2RAD
*wangle
);
1292 sprintf(buf
,"bonds that rotated more than %g degrees:\n"
1293 " atom 1 atom 2 angle previous, current, constraint length\n",
1295 fprintf(stderr
,"%s",buf
);
1298 fprintf(fplog
,"%s",buf
);
1301 for(b
=0;b
<ncons
;b
++)
1307 pbc_dx_aiuc(pbc
,x
[i
],x
[j
],v0
);
1308 pbc_dx_aiuc(pbc
,xprime
[i
],xprime
[j
],v1
);
1312 rvec_sub(x
[i
],x
[j
],v0
);
1313 rvec_sub(xprime
[i
],xprime
[j
],v1
);
1317 cosine
= iprod(v0
,v1
)/(d0
*d1
);
1320 sprintf(buf
," %6d %6d %5.1f %8.4f %8.4f %8.4f\n",
1321 ddglatnr(dd
,i
),ddglatnr(dd
,j
),
1322 RAD2DEG
*acos(cosine
),d0
,d1
,bllen
[b
]);
1323 fprintf(stderr
,"%s",buf
);
1326 fprintf(fplog
,"%s",buf
);
1328 if (!gmx_isfinite(d1
))
1330 gmx_fatal(FARGS
,"Bond length not finite.");
1336 if (*warncount
> maxwarn
)
1338 too_many_constraint_warnings(econtLINCS
,*warncount
);
1342 static void cconerr(gmx_domdec_t
*dd
,
1343 int ncons
,int *bla
,real
*bllen
,rvec
*x
,t_pbc
*pbc
,
1344 real
*ncons_loc
,real
*ssd
,real
*max
,int *imax
)
1346 real len
,d
,ma
,ssd2
,r2
;
1347 int *nlocat
,count
,b
,im
;
1350 if (dd
&& dd
->constraints
)
1352 nlocat
= dd_constraints_nlocalatoms(dd
);
1363 for(b
=0;b
<ncons
;b
++)
1367 pbc_dx_aiuc(pbc
,x
[bla
[2*b
]],x
[bla
[2*b
+1]],dx
);
1370 rvec_sub(x
[bla
[2*b
]],x
[bla
[2*b
+1]],dx
);
1373 len
= r2
*gmx_invsqrt(r2
);
1374 d
= fabs(len
/bllen
[b
]-1);
1375 if (d
> ma
&& (nlocat
==NULL
|| nlocat
[b
]))
1387 ssd2
+= nlocat
[b
]*d
*d
;
1392 *ncons_loc
= (nlocat
? 0.5 : 1)*count
;
1393 *ssd
= (nlocat
? 0.5 : 1)*ssd2
;
1398 static void dump_conf(gmx_domdec_t
*dd
,struct gmx_lincsdata
*li
,
1400 char *name
,gmx_bool bAll
,rvec
*x
,matrix box
)
1406 dd_get_constraint_range(dd
,&ac0
,&ac1
);
1408 sprintf(str
,"%s_%d_%d_%d.pdb",name
,dd
->ci
[XX
],dd
->ci
[YY
],dd
->ci
[ZZ
]);
1409 fp
= gmx_fio_fopen(str
,"w");
1410 fprintf(fp
,"CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1\n",
1411 10*norm(box
[XX
]),10*norm(box
[YY
]),10*norm(box
[ZZ
]),
1413 for(i
=0; i
<ac1
; i
++)
1415 if (i
< dd
->nat_home
|| (bAll
&& i
>= ac0
&& i
< ac1
))
1417 fprintf(fp
,"%-6s%5u %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
1418 "ATOM",ddglatnr(dd
,i
),"C","ALA",' ',i
+1,
1419 10*x
[i
][XX
],10*x
[i
][YY
],10*x
[i
][ZZ
],
1420 1.0,i
<dd
->nat_tot
? 0.0 : 1.0);
1425 for(i
=0; i
<li
->nc
; i
++)
1427 fprintf(fp
,"CONECT%5d%5d\n",
1428 ddglatnr(dd
,li
->bla
[2*i
]),
1429 ddglatnr(dd
,li
->bla
[2*i
+1]));
1435 gmx_bool
constrain_lincs(FILE *fplog
,gmx_bool bLog
,gmx_bool bEner
,
1437 gmx_large_int_t step
,
1438 struct gmx_lincsdata
*lincsd
,t_mdatoms
*md
,
1440 rvec
*x
,rvec
*xprime
,rvec
*min_proj
,
1441 matrix box
,t_pbc
*pbc
,
1442 real lambda
,real
*dvdlambda
,
1444 gmx_bool bCalcVir
,tensor rmdr
,
1447 int maxwarn
,int *warncount
)
1449 char buf
[STRLEN
],buf2
[22],buf3
[STRLEN
];
1450 int i
,warn
=0,p_imax
,error
;
1451 real ncons_loc
,p_ssd
,p_max
=0;
1457 if (lincsd
->nc
== 0 && cr
->dd
== NULL
)
1461 lincsd
->rmsd_data
[0] = 0;
1462 if (ir
->eI
== eiSD2
&& v
== NULL
)
1470 lincsd
->rmsd_data
[i
] = 0;
1476 if (econq
== econqCoord
)
1478 if (ir
->efep
!= efepNO
)
1480 if (md
->nMassPerturbed
&& lincsd
->matlam
!= md
->lambda
)
1482 set_lincs_matrix(lincsd
,md
->invmass
,md
->lambda
);
1485 for(i
=0; i
<lincsd
->nc
; i
++)
1487 lincsd
->bllen
[i
] = lincsd
->bllen0
[i
] + lambda
*lincsd
->ddist
[i
];
1491 if (lincsd
->ncg_flex
)
1493 /* Set the flexible constraint lengths to the old lengths */
1496 for(i
=0; i
<lincsd
->nc
; i
++)
1498 if (lincsd
->bllen
[i
] == 0) {
1499 pbc_dx_aiuc(pbc
,x
[lincsd
->bla
[2*i
]],x
[lincsd
->bla
[2*i
+1]],dx
);
1500 lincsd
->bllen
[i
] = norm(dx
);
1506 for(i
=0; i
<lincsd
->nc
; i
++)
1508 if (lincsd
->bllen
[i
] == 0)
1511 sqrt(distance2(x
[lincsd
->bla
[2*i
]],
1512 x
[lincsd
->bla
[2*i
+1]]));
1520 cconerr(cr
->dd
,lincsd
->nc
,lincsd
->bla
,lincsd
->bllen
,xprime
,pbc
,
1521 &ncons_loc
,&p_ssd
,&p_max
,&p_imax
);
1524 /* The (only) OpenMP parallel region of constrain_lincs */
1528 #pragma omp parallel for num_threads(lincsd->nth) schedule(static)
1529 for(th
=0; th
<lincsd
->nth
; th
++)
1531 do_lincs(x
,xprime
,box
,pbc
,lincsd
,th
,
1533 bCalcVir
|| (ir
->efep
!= efepNO
),
1534 ir
->LincsWarnAngle
,&warn
,
1535 invdt
,v
,bCalcVir
,rmdr
);
1539 if (ir
->efep
!= efepNO
)
1543 dt_2
= 1.0/(ir
->delta_t
*ir
->delta_t
);
1544 for(i
=0; (i
<lincsd
->nc
); i
++)
1546 dvdl
-= lincsd
->mlambda
[i
]*dt_2
*lincsd
->ddist
[i
];
1551 if (bLog
&& fplog
&& lincsd
->nc
> 0)
1553 fprintf(fplog
," Rel. Constraint Deviation: RMS MAX between atoms\n");
1554 fprintf(fplog
," Before LINCS %.6f %.6f %6d %6d\n",
1555 sqrt(p_ssd
/ncons_loc
),p_max
,
1556 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
]),
1557 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
+1]));
1561 cconerr(cr
->dd
,lincsd
->nc
,lincsd
->bla
,lincsd
->bllen
,xprime
,pbc
,
1562 &ncons_loc
,&p_ssd
,&p_max
,&p_imax
);
1563 /* Check if we are doing the second part of SD */
1564 if (ir
->eI
== eiSD2
&& v
== NULL
)
1572 lincsd
->rmsd_data
[0] = ncons_loc
;
1573 lincsd
->rmsd_data
[i
] = p_ssd
;
1577 lincsd
->rmsd_data
[0] = 0;
1578 lincsd
->rmsd_data
[1] = 0;
1579 lincsd
->rmsd_data
[2] = 0;
1581 if (bLog
&& fplog
&& lincsd
->nc
> 0)
1584 " After LINCS %.6f %.6f %6d %6d\n\n",
1585 sqrt(p_ssd
/ncons_loc
),p_max
,
1586 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
]),
1587 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
+1]));
1594 cconerr(cr
->dd
,lincsd
->nc
,lincsd
->bla
,lincsd
->bllen
,xprime
,pbc
,
1595 &ncons_loc
,&p_ssd
,&p_max
,&p_imax
);
1598 sprintf(buf3
," in simulation %d", cr
->ms
->sim
);
1604 sprintf(buf
,"\nStep %s, time %g (ps) LINCS WARNING%s\n"
1605 "relative constraint deviation after LINCS:\n"
1606 "rms %.6f, max %.6f (between atoms %d and %d)\n",
1607 gmx_step_str(step
,buf2
),ir
->init_t
+step
*ir
->delta_t
,
1609 sqrt(p_ssd
/ncons_loc
),p_max
,
1610 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
]),
1611 ddglatnr(cr
->dd
,lincsd
->bla
[2*p_imax
+1]));
1614 fprintf(fplog
,"%s",buf
);
1616 fprintf(stderr
,"%s",buf
);
1617 lincs_warning(fplog
,cr
->dd
,x
,xprime
,pbc
,
1618 lincsd
->nc
,lincsd
->bla
,lincsd
->bllen
,
1619 ir
->LincsWarnAngle
,maxwarn
,warncount
);
1621 bOK
= (p_max
< 0.5);
1624 if (lincsd
->ncg_flex
) {
1625 for(i
=0; (i
<lincsd
->nc
); i
++)
1626 if (lincsd
->bllen0
[i
] == 0 && lincsd
->ddist
[i
] == 0)
1627 lincsd
->bllen
[i
] = 0;
1632 do_lincsp(x
,xprime
,min_proj
,pbc
,lincsd
,md
->invmass
,econq
,dvdlambda
,
1636 /* count assuming nit=1 */
1637 inc_nrnb(nrnb
,eNR_LINCS
,lincsd
->nc
);
1638 inc_nrnb(nrnb
,eNR_LINCSMAT
,(2+lincsd
->nOrder
)*lincsd
->ncc
);
1639 if (lincsd
->ntriangle
> 0)
1641 inc_nrnb(nrnb
,eNR_LINCSMAT
,lincsd
->nOrder
*lincsd
->ncc_triangle
);
1645 inc_nrnb(nrnb
,eNR_CONSTR_V
,lincsd
->nc
*2);
1649 inc_nrnb(nrnb
,eNR_CONSTR_VIR
,lincsd
->nc
);