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
53 #include "mtop_util.h"
55 /* Routines to send/recieve coordinates and force
56 * of constructing atoms.
59 static void move_construct_x(t_comm_vsites *vsitecomm, rvec x[], t_commrec *cr)
65 sendbuf = vsitecomm->send_buf;
66 recvbuf = vsitecomm->recv_buf;
69 /* Prepare pulse left by copying to send buffer */
70 for(i=0;i<vsitecomm->left_export_nconstruct;i++)
72 ia = vsitecomm->left_export_construct[i];
73 copy_rvec(x[ia],sendbuf[i]);
76 /* Pulse coordinates left */
77 gmx_tx_rx_real(cr,GMX_LEFT,(real *)sendbuf,3*vsitecomm->left_export_nconstruct,GMX_RIGHT,(real *)recvbuf,3*vsitecomm->right_import_nconstruct);
79 /* Copy from receive buffer to coordinate array */
80 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
82 ia = vsitecomm->right_import_construct[i];
83 copy_rvec(recvbuf[i],x[ia]);
86 /* Prepare pulse right by copying to send buffer */
87 for(i=0;i<vsitecomm->right_export_nconstruct;i++)
89 ia = vsitecomm->right_export_construct[i];
90 copy_rvec(x[ia],sendbuf[i]);
93 /* Pulse coordinates right */
94 gmx_tx_rx_real(cr,GMX_RIGHT,(real *)sendbuf,3*vsitecomm->right_export_nconstruct,GMX_LEFT,(real *)recvbuf,3*vsitecomm->left_import_nconstruct);
96 /* Copy from receive buffer to coordinate array */
97 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
99 ia = vsitecomm->left_import_construct[i];
100 copy_rvec(recvbuf[i],x[ia]);
105 static void move_construct_f(t_comm_vsites *vsitecomm, rvec f[], t_commrec *cr)
111 sendbuf = vsitecomm->send_buf;
112 recvbuf = vsitecomm->recv_buf;
114 /* Prepare pulse right by copying to send buffer */
115 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
117 ia = vsitecomm->right_import_construct[i];
118 copy_rvec(f[ia],sendbuf[i]);
119 clear_rvec(f[ia]); /* Zero it here after moving, just to simplify debug book-keeping... */
122 /* Pulse forces right */
123 gmx_tx_rx_real(cr,GMX_RIGHT,(real *)sendbuf,3*vsitecomm->right_import_nconstruct,GMX_LEFT,(real *)recvbuf,3*vsitecomm->left_export_nconstruct);
125 /* Copy from receive buffer to coordinate array */
126 for(i=0;i<vsitecomm->left_export_nconstruct;i++)
128 ia = vsitecomm->left_export_construct[i];
129 rvec_inc(f[ia],recvbuf[i]);
132 /* Prepare pulse left by copying to send buffer */
133 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
135 ia = vsitecomm->left_import_construct[i];
136 copy_rvec(f[ia],sendbuf[i]);
137 clear_rvec(f[ia]); /* Zero it here after moving, just to simplify debug book-keeping... */
140 /* Pulse coordinates left */
141 gmx_tx_rx_real(cr,GMX_LEFT,(real *)sendbuf,3*vsitecomm->left_import_nconstruct,GMX_RIGHT,(real *)recvbuf,3*vsitecomm->right_export_nconstruct);
143 /* Copy from receive buffer to coordinate array */
144 for(i=0;i<vsitecomm->right_export_nconstruct;i++)
146 ia = vsitecomm->right_export_construct[i];
147 rvec_inc(f[ia],recvbuf[i]);
150 /* All forces are now on the home processors */
155 pd_clear_nonlocal_constructs(t_comm_vsites *vsitecomm, rvec f[])
159 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
161 ia = vsitecomm->left_import_construct[i];
164 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
166 ia = vsitecomm->right_import_construct[i];
173 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
176 return pbc_dx_aiuc(pbc,xi,xj,dx);
184 /* Vsite construction routines */
186 static void constr_vsite2(rvec xi,rvec xj,rvec x,real a,t_pbc *pbc)
195 pbc_dx_aiuc(pbc,xj,xi,dx);
196 x[XX] = xi[XX] + a*dx[XX];
197 x[YY] = xi[YY] + a*dx[YY];
198 x[ZZ] = xi[ZZ] + a*dx[ZZ];
200 x[XX] = b*xi[XX] + a*xj[XX];
201 x[YY] = b*xi[YY] + a*xj[YY];
202 x[ZZ] = b*xi[ZZ] + a*xj[ZZ];
206 /* TOTAL: 10 flops */
209 static void constr_vsite3(rvec xi,rvec xj,rvec xk,rvec x,real a,real b,
219 pbc_dx_aiuc(pbc,xj,xi,dxj);
220 pbc_dx_aiuc(pbc,xk,xi,dxk);
221 x[XX] = xi[XX] + a*dxj[XX] + b*dxk[XX];
222 x[YY] = xi[YY] + a*dxj[YY] + b*dxk[YY];
223 x[ZZ] = xi[ZZ] + a*dxj[ZZ] + b*dxk[ZZ];
225 x[XX] = c*xi[XX] + a*xj[XX] + b*xk[XX];
226 x[YY] = c*xi[YY] + a*xj[YY] + b*xk[YY];
227 x[ZZ] = c*xi[ZZ] + a*xj[ZZ] + b*xk[ZZ];
231 /* TOTAL: 17 flops */
234 static void constr_vsite3FD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b,
240 pbc_rvec_sub(pbc,xj,xi,xij);
241 pbc_rvec_sub(pbc,xk,xj,xjk);
244 /* temp goes from i to a point on the line jk */
245 temp[XX] = xij[XX] + a*xjk[XX];
246 temp[YY] = xij[YY] + a*xjk[YY];
247 temp[ZZ] = xij[ZZ] + a*xjk[ZZ];
250 c=b*gmx_invsqrt(iprod(temp,temp));
253 x[XX] = xi[XX] + c*temp[XX];
254 x[YY] = xi[YY] + c*temp[YY];
255 x[ZZ] = xi[ZZ] + c*temp[ZZ];
258 /* TOTAL: 34 flops */
261 static void constr_vsite3FAD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b, t_pbc *pbc)
264 real a1,b1,c1,invdij;
266 pbc_rvec_sub(pbc,xj,xi,xij);
267 pbc_rvec_sub(pbc,xk,xj,xjk);
270 invdij = gmx_invsqrt(iprod(xij,xij));
271 c1 = invdij * invdij * iprod(xij,xjk);
272 xp[XX] = xjk[XX] - c1*xij[XX];
273 xp[YY] = xjk[YY] - c1*xij[YY];
274 xp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
276 b1 = b*gmx_invsqrt(iprod(xp,xp));
279 x[XX] = xi[XX] + a1*xij[XX] + b1*xp[XX];
280 x[YY] = xi[YY] + a1*xij[YY] + b1*xp[YY];
281 x[ZZ] = xi[ZZ] + a1*xij[ZZ] + b1*xp[ZZ];
284 /* TOTAL: 63 flops */
287 static void constr_vsite3OUT(rvec xi,rvec xj,rvec xk,rvec x,
288 real a,real b,real c,t_pbc *pbc)
292 pbc_rvec_sub(pbc,xj,xi,xij);
293 pbc_rvec_sub(pbc,xk,xi,xik);
297 x[XX] = xi[XX] + a*xij[XX] + b*xik[XX] + c*temp[XX];
298 x[YY] = xi[YY] + a*xij[YY] + b*xik[YY] + c*temp[YY];
299 x[ZZ] = xi[ZZ] + a*xij[ZZ] + b*xik[ZZ] + c*temp[ZZ];
302 /* TOTAL: 33 flops */
305 static void constr_vsite4FD(rvec xi,rvec xj,rvec xk,rvec xl,rvec x,
306 real a,real b,real c,t_pbc *pbc)
308 rvec xij,xjk,xjl,temp;
311 pbc_rvec_sub(pbc,xj,xi,xij);
312 pbc_rvec_sub(pbc,xk,xj,xjk);
313 pbc_rvec_sub(pbc,xl,xj,xjl);
316 /* temp goes from i to a point on the plane jkl */
317 temp[XX] = xij[XX] + a*xjk[XX] + b*xjl[XX];
318 temp[YY] = xij[YY] + a*xjk[YY] + b*xjl[YY];
319 temp[ZZ] = xij[ZZ] + a*xjk[ZZ] + b*xjl[ZZ];
322 d=c*gmx_invsqrt(iprod(temp,temp));
325 x[XX] = xi[XX] + d*temp[XX];
326 x[YY] = xi[YY] + d*temp[YY];
327 x[ZZ] = xi[ZZ] + d*temp[ZZ];
330 /* TOTAL: 43 flops */
334 static void constr_vsite4FDN(rvec xi,rvec xj,rvec xk,rvec xl,rvec x,
335 real a,real b,real c,t_pbc *pbc)
337 rvec xij,xik,xil,ra,rb,rja,rjb,rm;
340 pbc_rvec_sub(pbc,xj,xi,xij);
341 pbc_rvec_sub(pbc,xk,xi,xik);
342 pbc_rvec_sub(pbc,xl,xi,xil);
355 rvec_sub(ra,xij,rja);
356 rvec_sub(rb,xij,rjb);
362 d=c*gmx_invsqrt(norm2(rm));
365 x[XX] = xi[XX] + d*rm[XX];
366 x[YY] = xi[YY] + d*rm[YY];
367 x[ZZ] = xi[ZZ] + d*rm[ZZ];
370 /* TOTAL: 47 flops */
374 static int constr_vsiten(t_iatom *ia, t_iparams ip[],
382 n3 = 3*ip[ia[0]].vsiten.n;
387 for(i=3; i<n3; i+=3) {
389 a = ip[ia[i]].vsiten.a;
391 pbc_dx_aiuc(pbc,x[ai],x1,dx);
393 rvec_sub(x[ai],x1,dx);
395 dsum[XX] += a*dx[XX];
396 dsum[YY] += a*dx[YY];
397 dsum[ZZ] += a*dx[ZZ];
401 x[av][XX] = x1[XX] + dsum[XX];
402 x[av][YY] = x1[YY] + dsum[YY];
403 x[av][ZZ] = x1[ZZ] + dsum[ZZ];
409 void construct_vsites(FILE *log,gmx_vsite_t *vsite,
410 rvec x[],t_nrnb *nrnb,
412 t_iparams ip[],t_ilist ilist[],
413 int ePBC,gmx_bool bMolPBC,t_graph *graph,
414 t_commrec *cr,matrix box)
417 real a1,b1,c1,inv_dt;
418 int i,inc,ii,nra,nr,tp,ftype;
419 t_iatom avsite,ai,aj,ak,al,pbc_atom;
421 t_pbc pbc,*pbc_null,*pbc_null2;
423 int *vsite_pbc,ishift;
424 rvec reftmp,vtmp,rtmp;
426 bDomDec = cr && DOMAINDECOMP(cr);
428 /* We only need to do pbc when we have inter-cg vsites */
429 if (ePBC != epbcNONE && (bDomDec || bMolPBC) && vsite->n_intercg_vsite) {
430 /* This is wasting some CPU time as we now do this multiple times
431 * per MD step. But how often do we have vsites with full pbc?
433 pbc_null = set_pbc_dd(&pbc,ePBC,cr!=NULL ? cr->dd : NULL,FALSE,box);
440 dd_move_x_vsites(cr->dd,box,x);
441 } else if (vsite->bPDvsitecomm) {
442 /* I'm not sure whether the periodicity and shift are guaranteed
443 * to be consistent between different nodes when running e.g. polymers
444 * in parallel. In this special case we thus unshift/shift,
445 * but only when necessary. This is to make sure the coordinates
446 * we move don't end up a box away...
449 unshift_self(graph,box,x);
451 move_construct_x(vsite->vsitecomm,x,cr);
454 shift_self(graph,box,x);
465 for(ftype=0; (ftype<F_NRE); ftype++) {
466 if (interaction_function[ftype].flags & IF_VSITE) {
467 nra = interaction_function[ftype].nratoms;
468 nr = ilist[ftype].nr;
469 ia = ilist[ftype].iatoms;
472 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
480 if (ftype != idef->functype[tp])
481 gmx_incons("Function types for vsites wrong");
484 /* The vsite and constructing atoms */
489 /* Constants for constructing vsites */
491 /* Check what kind of pbc we need to use */
493 pbc_atom = vsite_pbc[i/(1+nra)];
496 /* We need to copy the coordinates here,
497 * single for single atom cg's pbc_atom is the vsite itself.
499 copy_rvec(x[pbc_atom],xpbc);
501 pbc_null2 = pbc_null;
508 /* Copy the old position */
509 copy_rvec(x[avsite],xv);
511 /* Construct the vsite depending on type */
515 constr_vsite2(x[ai],x[aj],x[avsite],a1,pbc_null2);
520 constr_vsite3(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
525 constr_vsite3FD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
530 constr_vsite3FAD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
536 constr_vsite3OUT(x[ai],x[aj],x[ak],x[avsite],a1,b1,c1,pbc_null2);
543 constr_vsite4FD(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1,
551 constr_vsite4FDN(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1,
555 inc = constr_vsiten(ia,ip,x,pbc_null2);
558 gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
559 ftype,__FILE__,__LINE__);
563 /* Match the pbc of this vsite to the rest of its charge group */
564 ishift = pbc_dx_aiuc(pbc_null,x[avsite],xpbc,dx);
565 if (ishift != CENTRAL)
566 rvec_add(xpbc,dx,x[avsite]);
569 /* Calculate velocity of vsite... */
570 rvec_sub(x[avsite],xv,vv);
571 svmul(inv_dt,vv,v[avsite]);
574 /* Increment loop variables */
582 static void spread_vsite2(t_iatom ia[],real a,
583 rvec x[],rvec f[],rvec fshift[],
584 t_pbc *pbc,t_graph *g)
605 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,av),di);
607 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di);
610 siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx);
611 sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
617 if (fshift && (siv != CENTRAL || sij != CENTRAL)) {
618 rvec_inc(fshift[siv],f[av]);
619 rvec_dec(fshift[CENTRAL],fi);
620 rvec_dec(fshift[sij],fj);
623 /* TOTAL: 13 flops */
626 void construct_vsites_mtop(FILE *log,gmx_vsite_t *vsite,
627 gmx_mtop_t *mtop,rvec x[])
630 gmx_molblock_t *molb;
634 for(mb=0; mb<mtop->nmolblock; mb++) {
635 molb = &mtop->molblock[mb];
636 molt = &mtop->moltype[molb->type];
637 for(mol=0; mol<molb->nmol; mol++) {
638 construct_vsites(log,vsite,x+as,NULL,0.0,NULL,
639 mtop->ffparams.iparams,molt->ilist,
640 epbcNONE,TRUE,NULL,NULL,NULL);
641 as += molt->atoms.nr;
646 static void spread_vsite3(t_iatom ia[],real a,real b,
647 rvec x[],rvec f[],rvec fshift[],
648 t_pbc *pbc,t_graph *g)
660 svmul(1-a-b,f[av],fi);
671 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ia[1]),di);
673 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di);
675 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),di);
678 siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx);
679 sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
680 sik = pbc_dx_aiuc(pbc,x[ai],x[ak],dx);
687 if (fshift && (siv!=CENTRAL || sij!=CENTRAL || sik!=CENTRAL)) {
688 rvec_inc(fshift[siv],f[av]);
689 rvec_dec(fshift[CENTRAL],fi);
690 rvec_dec(fshift[sij],fj);
691 rvec_dec(fshift[sik],fk);
694 /* TOTAL: 20 flops */
697 static void spread_vsite3FD(t_iatom ia[],real a,real b,
698 rvec x[],rvec f[],rvec fshift[],
699 gmx_bool VirCorr,matrix dxdf,
700 t_pbc *pbc,t_graph *g)
702 real fx,fy,fz,c,invl,fproj,a1;
703 rvec xvi,xij,xjk,xix,fv,temp;
714 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
715 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
718 /* xix goes from i to point x on the line jk */
719 xix[XX]=xij[XX]+a*xjk[XX];
720 xix[YY]=xij[YY]+a*xjk[YY];
721 xix[ZZ]=xij[ZZ]+a*xjk[ZZ];
724 invl=gmx_invsqrt(iprod(xix,xix));
728 fproj=iprod(xix,fv)*invl*invl; /* = (xix . f)/(xix . xix) */
730 temp[XX]=c*(fv[XX]-fproj*xix[XX]);
731 temp[YY]=c*(fv[YY]-fproj*xix[YY]);
732 temp[ZZ]=c*(fv[ZZ]-fproj*xix[ZZ]);
735 /* c is already calculated in constr_vsite3FD
736 storing c somewhere will save 26 flops! */
739 f[ai][XX] += fv[XX] - temp[XX];
740 f[ai][YY] += fv[YY] - temp[YY];
741 f[ai][ZZ] += fv[ZZ] - temp[ZZ];
742 f[aj][XX] += a1*temp[XX];
743 f[aj][YY] += a1*temp[YY];
744 f[aj][ZZ] += a1*temp[ZZ];
745 f[ak][XX] += a*temp[XX];
746 f[ak][YY] += a*temp[YY];
747 f[ak][ZZ] += a*temp[ZZ];
751 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
753 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
755 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
758 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
763 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL)) {
764 rvec_dec(fshift[svi],fv);
765 fshift[CENTRAL][XX] += fv[XX] - (1 + a)*temp[XX];
766 fshift[CENTRAL][YY] += fv[YY] - (1 + a)*temp[YY];
767 fshift[CENTRAL][ZZ] += fv[ZZ] - (1 + a)*temp[ZZ];
768 fshift[ sji][XX] += temp[XX];
769 fshift[ sji][YY] += temp[YY];
770 fshift[ sji][ZZ] += temp[ZZ];
771 fshift[ skj][XX] += a*temp[XX];
772 fshift[ skj][YY] += a*temp[YY];
773 fshift[ skj][ZZ] += a*temp[ZZ];
778 /* When VirCorr=TRUE, the virial for the current forces is not
779 * calculated from the redistributed forces. This means that
780 * the effect of non-linear virtual site constructions on the virial
781 * needs to be added separately. This contribution can be calculated
782 * in many ways, but the simplest and cheapest way is to use
783 * the first constructing atom ai as a reference position in space:
784 * subtract (xv-xi)*fv and add (xj-xi)*fj + (xk-xi)*fk.
789 pbc_rvec_sub(pbc,x[av],x[ai],xiv);
795 /* As xix is a linear combination of j and k, use that here */
796 dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
801 /* TOTAL: 61 flops */
804 static void spread_vsite3FAD(t_iatom ia[],real a,real b,
805 rvec x[],rvec f[],rvec fshift[],
806 gmx_bool VirCorr,matrix dxdf,
807 t_pbc *pbc,t_graph *g)
809 rvec xvi,xij,xjk,xperp,Fpij,Fppp,fv,f1,f2,f3;
810 real a1,b1,c1,c2,invdij,invdij2,invdp,fproj;
819 copy_rvec(f[ia[1]],fv);
821 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
822 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
825 invdij = gmx_invsqrt(iprod(xij,xij));
826 invdij2 = invdij * invdij;
827 c1 = iprod(xij,xjk) * invdij2;
828 xperp[XX] = xjk[XX] - c1*xij[XX];
829 xperp[YY] = xjk[YY] - c1*xij[YY];
830 xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
831 /* xperp in plane ijk, perp. to ij */
832 invdp = gmx_invsqrt(iprod(xperp,xperp));
837 /* a1, b1 and c1 are already calculated in constr_vsite3FAD
838 storing them somewhere will save 45 flops! */
840 fproj=iprod(xij ,fv)*invdij2;
841 svmul(fproj, xij, Fpij); /* proj. f on xij */
842 svmul(iprod(xperp,fv)*invdp*invdp,xperp,Fppp); /* proj. f on xperp */
843 svmul(b1*fproj, xperp,f3);
846 rvec_sub(fv,Fpij,f1); /* f1 = f - Fpij */
847 rvec_sub(f1,Fppp,f2); /* f2 = f - Fpij - Fppp */
848 for (d=0; (d<DIM); d++) {
855 f[ai][XX] += fv[XX] - f1[XX] + c1*f2[XX] + f3[XX];
856 f[ai][YY] += fv[YY] - f1[YY] + c1*f2[YY] + f3[YY];
857 f[ai][ZZ] += fv[ZZ] - f1[ZZ] + c1*f2[ZZ] + f3[ZZ];
858 f[aj][XX] += f1[XX] - c2*f2[XX] - f3[XX];
859 f[aj][YY] += f1[YY] - c2*f2[YY] - f3[YY];
860 f[aj][ZZ] += f1[ZZ] - c2*f2[ZZ] - f3[ZZ];
867 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
869 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
871 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
874 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
879 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL)) {
880 rvec_dec(fshift[svi],fv);
881 fshift[CENTRAL][XX] += fv[XX] - f1[XX] - (1-c1)*f2[XX] + f3[XX];
882 fshift[CENTRAL][YY] += fv[YY] - f1[YY] - (1-c1)*f2[YY] + f3[YY];
883 fshift[CENTRAL][ZZ] += fv[ZZ] - f1[ZZ] - (1-c1)*f2[ZZ] + f3[ZZ];
884 fshift[ sji][XX] += f1[XX] - c1 *f2[XX] - f3[XX];
885 fshift[ sji][YY] += f1[YY] - c1 *f2[YY] - f3[YY];
886 fshift[ sji][ZZ] += f1[ZZ] - c1 *f2[ZZ] - f3[ZZ];
887 fshift[ skj][XX] += f2[XX];
888 fshift[ skj][YY] += f2[YY];
889 fshift[ skj][ZZ] += f2[ZZ];
897 pbc_rvec_sub(pbc,x[av],x[ai],xiv);
903 /* Note that xik=xij+xjk, so we have to add xij*f2 */
906 + xij[i]*(f1[j] + (1 - c2)*f2[j] - f3[j])
912 /* TOTAL: 113 flops */
915 static void spread_vsite3OUT(t_iatom ia[],real a,real b,real c,
916 rvec x[],rvec f[],rvec fshift[],
917 gmx_bool VirCorr,matrix dxdf,
918 t_pbc *pbc,t_graph *g)
920 rvec xvi,xij,xik,fv,fj,fk;
931 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
932 ski = pbc_rvec_sub(pbc,x[ak],x[ai],xik);
942 fj[XX] = a*fv[XX] - xik[ZZ]*cfy + xik[YY]*cfz;
943 fj[YY] = xik[ZZ]*cfx + a*fv[YY] - xik[XX]*cfz;
944 fj[ZZ] = -xik[YY]*cfx + xik[XX]*cfy + a*fv[ZZ];
946 fk[XX] = b*fv[XX] + xij[ZZ]*cfy - xij[YY]*cfz;
947 fk[YY] = -xij[ZZ]*cfx + b*fv[YY] + xij[XX]*cfz;
948 fk[ZZ] = xij[YY]*cfx - xij[XX]*cfy + b*fv[ZZ];
951 f[ai][XX] += fv[XX] - fj[XX] - fk[XX];
952 f[ai][YY] += fv[YY] - fj[YY] - fk[YY];
953 f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
959 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
961 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
963 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,ai),di);
966 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
971 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || ski!=CENTRAL)) {
972 rvec_dec(fshift[svi],fv);
973 fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX];
974 fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY];
975 fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
976 rvec_inc(fshift[sji],fj);
977 rvec_inc(fshift[ski],fk);
985 pbc_rvec_sub(pbc,x[av],x[ai],xiv);
991 dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j];
996 /* TOTAL: 54 flops */
999 static void spread_vsite4FD(t_iatom ia[],real a,real b,real c,
1000 rvec x[],rvec f[],rvec fshift[],
1001 gmx_bool VirCorr,matrix dxdf,
1002 t_pbc *pbc,t_graph *g)
1004 real d,invl,fproj,a1;
1005 rvec xvi,xij,xjk,xjl,xix,fv,temp;
1006 atom_id av,ai,aj,ak,al;
1008 int svi,sji,skj,slj,m;
1016 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
1017 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
1018 slj = pbc_rvec_sub(pbc,x[al],x[aj],xjl);
1021 /* xix goes from i to point x on the plane jkl */
1022 for(m=0; m<DIM; m++)
1023 xix[m] = xij[m] + a*xjk[m] + b*xjl[m];
1026 invl=gmx_invsqrt(iprod(xix,xix));
1028 /* 4 + ?10? flops */
1030 copy_rvec(f[av],fv);
1032 fproj=iprod(xix,fv)*invl*invl; /* = (xix . f)/(xix . xix) */
1034 for(m=0; m<DIM; m++)
1035 temp[m] = d*(fv[m] - fproj*xix[m]);
1038 /* c is already calculated in constr_vsite3FD
1039 storing c somewhere will save 35 flops! */
1042 for(m=0; m<DIM; m++) {
1043 f[ai][m] += fv[m] - temp[m];
1044 f[aj][m] += a1*temp[m];
1045 f[ak][m] += a*temp[m];
1046 f[al][m] += b*temp[m];
1051 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
1053 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
1055 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
1057 ivec_sub(SHIFT_IVEC(g,al),SHIFT_IVEC(g,aj),di);
1060 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
1066 (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL || slj!=CENTRAL)) {
1067 rvec_dec(fshift[svi],fv);
1068 for(m=0; m<DIM; m++) {
1069 fshift[CENTRAL][m] += fv[m] - (1 + a + b)*temp[m];
1070 fshift[ sji][m] += temp[m];
1071 fshift[ skj][m] += a*temp[m];
1072 fshift[ slj][m] += b*temp[m];
1081 pbc_rvec_sub(pbc,x[av],x[ai],xiv);
1083 for(i=0; i<DIM; i++)
1085 for(j=0; j<DIM; j++)
1087 dxdf[i][j] += -xiv[i]*fv[j] + xix[i]*temp[j];
1092 /* TOTAL: 77 flops */
1096 static void spread_vsite4FDN(t_iatom ia[],real a,real b,real c,
1097 rvec x[],rvec f[],rvec fshift[],
1098 gmx_bool VirCorr,matrix dxdf,
1099 t_pbc *pbc,t_graph *g)
1101 rvec xvi,xij,xik,xil,ra,rb,rja,rjb,rab,rm,rt;
1107 int svi,sij,sik,sil;
1109 /* DEBUG: check atom indices */
1116 copy_rvec(f[av],fv);
1118 sij = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
1119 sik = pbc_rvec_sub(pbc,x[ak],x[ai],xik);
1120 sil = pbc_rvec_sub(pbc,x[al],x[ai],xil);
1133 rvec_sub(ra,xij,rja);
1134 rvec_sub(rb,xij,rjb);
1135 rvec_sub(rb,ra,rab);
1141 invrm=gmx_invsqrt(norm2(rm));
1145 cfx = c*invrm*fv[XX];
1146 cfy = c*invrm*fv[YY];
1147 cfz = c*invrm*fv[ZZ];
1158 fj[XX] = ( -rm[XX]*rt[XX]) * cfx + ( rab[ZZ]-rm[YY]*rt[XX]) * cfy + (-rab[YY]-rm[ZZ]*rt[XX]) * cfz;
1159 fj[YY] = (-rab[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + ( rab[XX]-rm[ZZ]*rt[YY]) * cfz;
1160 fj[ZZ] = ( rab[YY]-rm[XX]*rt[ZZ]) * cfx + (-rab[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1171 fk[XX] = ( -rm[XX]*rt[XX]) * cfx + (-a*rjb[ZZ]-rm[YY]*rt[XX]) * cfy + ( a*rjb[YY]-rm[ZZ]*rt[XX]) * cfz;
1172 fk[YY] = ( a*rjb[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + (-a*rjb[XX]-rm[ZZ]*rt[YY]) * cfz;
1173 fk[ZZ] = (-a*rjb[YY]-rm[XX]*rt[ZZ]) * cfx + ( a*rjb[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1184 fl[XX] = ( -rm[XX]*rt[XX]) * cfx + ( b*rja[ZZ]-rm[YY]*rt[XX]) * cfy + (-b*rja[YY]-rm[ZZ]*rt[XX]) * cfz;
1185 fl[YY] = (-b*rja[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + ( b*rja[XX]-rm[ZZ]*rt[YY]) * cfz;
1186 fl[ZZ] = ( b*rja[YY]-rm[XX]*rt[ZZ]) * cfx + (-b*rja[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1189 f[ai][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1190 f[ai][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1191 f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1198 ivec_sub(SHIFT_IVEC(g,av),SHIFT_IVEC(g,ai),di);
1200 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
1202 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,ai),di);
1204 ivec_sub(SHIFT_IVEC(g,al),SHIFT_IVEC(g,ai),di);
1207 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
1212 if (fshift && (svi!=CENTRAL || sij!=CENTRAL || sik!=CENTRAL || sil!=CENTRAL)) {
1213 rvec_dec(fshift[svi],fv);
1214 fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1215 fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1216 fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1217 rvec_inc(fshift[sij],fj);
1218 rvec_inc(fshift[sik],fk);
1219 rvec_inc(fshift[sil],fl);
1227 pbc_rvec_sub(pbc,x[av],x[ai],xiv);
1229 for(i=0; i<DIM; i++)
1231 for(j=0; j<DIM; j++)
1233 dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j] + xik[i]*fk[j] + xil[i]*fl[j];
1238 /* Total: 207 flops (Yuck!) */
1242 static int spread_vsiten(t_iatom ia[],t_iparams ip[],
1243 rvec x[],rvec f[],rvec fshift[],
1244 t_pbc *pbc,t_graph *g)
1252 n3 = 3*ip[ia[0]].vsiten.n;
1254 copy_rvec(x[av],xv);
1256 for(i=0; i<n3; i+=3) {
1259 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,av),di);
1262 siv = pbc_dx_aiuc(pbc,x[ai],xv,dx);
1266 a = ip[ia[i]].vsiten.a;
1269 if (fshift && siv != CENTRAL) {
1270 rvec_inc(fshift[siv],fi);
1271 rvec_dec(fshift[CENTRAL],fi);
1280 void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
1281 rvec x[],rvec f[],rvec *fshift,
1282 gmx_bool VirCorr,matrix vir,
1283 t_nrnb *nrnb,t_idef *idef,
1284 int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
1288 int i,inc,m,nra,nr,tp,ftype;
1289 int nd2,nd3,nd3FD,nd3FAD,nd3OUT,nd4FD,nd4FDN,ndN;
1292 t_pbc pbc,*pbc_null,*pbc_null2;
1296 /* We only need to do pbc when we have inter-cg vsites */
1297 if ((DOMAINDECOMP(cr) || bMolPBC) && vsite->n_intercg_vsite) {
1298 /* This is wasting some CPU time as we now do this multiple times
1299 * per MD step. But how often do we have vsites with full pbc?
1301 pbc_null = set_pbc_dd(&pbc,ePBC,cr->dd,FALSE,box);
1306 if (DOMAINDECOMP(cr))
1308 dd_clear_f_vsites(cr->dd,f);
1310 else if (PARTDECOMP(cr) && vsite->vsitecomm != NULL)
1312 pd_clear_nonlocal_constructs(vsite->vsitecomm,f);
1331 /* this loop goes backwards to be able to build *
1332 * higher type vsites from lower types */
1334 for(ftype=F_NRE-1; (ftype>=0); ftype--) {
1335 if (interaction_function[ftype].flags & IF_VSITE) {
1336 nra = interaction_function[ftype].nratoms;
1337 nr = idef->il[ftype].nr;
1338 ia = idef->il[ftype].iatoms;
1341 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
1346 for(i=0; (i<nr); ) {
1347 /* Check if we need to apply pbc for this vsite */
1349 if (vsite_pbc[i/(1+nra)] > -2)
1350 pbc_null2 = pbc_null;
1356 if (ftype != idef->functype[tp])
1357 gmx_incons("Functiontypes for vsites wrong");
1359 /* Constants for constructing */
1360 a1 = ip[tp].vsite.a;
1361 /* Construct the vsite depending on type */
1365 spread_vsite2(ia,a1,x,f,fshift,pbc_null2,g);
1369 b1 = ip[tp].vsite.b;
1370 spread_vsite3(ia,a1,b1,x,f,fshift,pbc_null2,g);
1374 b1 = ip[tp].vsite.b;
1375 spread_vsite3FD(ia,a1,b1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
1379 b1 = ip[tp].vsite.b;
1380 spread_vsite3FAD(ia,a1,b1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
1384 b1 = ip[tp].vsite.b;
1385 c1 = ip[tp].vsite.c;
1386 spread_vsite3OUT(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
1390 b1 = ip[tp].vsite.b;
1391 c1 = ip[tp].vsite.c;
1392 spread_vsite4FD(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
1396 b1 = ip[tp].vsite.b;
1397 c1 = ip[tp].vsite.c;
1398 spread_vsite4FDN(ia,a1,b1,c1,x,f,fshift,VirCorr,dxdf,pbc_null2,g);
1402 inc = spread_vsiten(ia,ip,x,f,fshift,pbc_null2,g);
1406 gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
1407 ftype,__FILE__,__LINE__);
1409 clear_rvec(f[ia[1]]);
1411 /* Increment loop variables */
1422 for(i=0; i<DIM; i++)
1424 for(j=0; j<DIM; j++)
1426 vir[i][j] += -0.5*dxdf[i][j];
1431 inc_nrnb(nrnb,eNR_VSITE2, nd2 );
1432 inc_nrnb(nrnb,eNR_VSITE3, nd3 );
1433 inc_nrnb(nrnb,eNR_VSITE3FD, nd3FD );
1434 inc_nrnb(nrnb,eNR_VSITE3FAD,nd3FAD );
1435 inc_nrnb(nrnb,eNR_VSITE3OUT,nd3OUT );
1436 inc_nrnb(nrnb,eNR_VSITE4FD, nd4FD );
1437 inc_nrnb(nrnb,eNR_VSITE4FDN,nd4FDN );
1438 inc_nrnb(nrnb,eNR_VSITEN, ndN );
1440 if (DOMAINDECOMP(cr)) {
1441 dd_move_f_vsites(cr->dd,f,fshift);
1442 } else if (vsite->bPDvsitecomm) {
1443 /* We only move forces here, and they are independent of shifts */
1444 move_construct_f(vsite->vsitecomm,f,cr);
1448 static int *atom2cg(t_block *cgs)
1452 snew(a2cg,cgs->index[cgs->nr]);
1453 for(cg=0; cg<cgs->nr; cg++) {
1454 for(i=cgs->index[cg]; i<cgs->index[cg+1]; i++)
1461 static int count_intercg_vsite(gmx_mtop_t *mtop)
1463 int mb,mt,ftype,nral,i,cg,a;
1464 gmx_molblock_t *molb;
1465 gmx_moltype_t *molt;
1469 int n_intercg_vsite;
1471 n_intercg_vsite = 0;
1472 for(mb=0; mb<mtop->nmolblock; mb++) {
1473 molb = &mtop->molblock[mb];
1474 molt = &mtop->moltype[molb->type];
1475 a2cg = atom2cg(&molt->cgs);
1476 for(ftype=0; ftype<F_NRE; ftype++) {
1477 if (interaction_function[ftype].flags & IF_VSITE) {
1479 il = &molt->ilist[ftype];
1481 for(i=0; i<il->nr; i+=1+nral) {
1483 for(a=1; a<nral; a++) {
1484 if (a2cg[ia[1+a]] != cg) {
1485 n_intercg_vsite += molb->nmol;
1495 return n_intercg_vsite;
1498 static int **get_vsite_pbc(t_iparams *iparams,t_ilist *ilist,
1499 t_atom *atom,t_mdatoms *md,
1500 t_block *cgs,int *a2cg)
1502 int ftype,nral,i,j,vsi,vsite,cg_v,cg_c,a,nc3=0;
1505 int **vsite_pbc,*vsite_pbc_f;
1507 gmx_bool bViteOnlyCG_and_FirstAtom;
1509 /* Make an array that tells if the pbc of an atom is set */
1510 snew(pbc_set,cgs->index[cgs->nr]);
1511 /* PBC is set for all non vsites */
1512 for(a=0; a<cgs->index[cgs->nr]; a++) {
1513 if ((atom && atom[a].ptype != eptVSite) ||
1514 (md && md->ptype[a] != eptVSite)) {
1519 snew(vsite_pbc,F_VSITEN-F_VSITE2+1);
1521 for(ftype=0; ftype<F_NRE; ftype++) {
1522 if (interaction_function[ftype].flags & IF_VSITE) {
1527 snew(vsite_pbc[ftype-F_VSITE2],il->nr/(1+nral));
1528 vsite_pbc_f = vsite_pbc[ftype-F_VSITE2];
1531 while (i < il->nr) {
1535 /* A value of -2 signals that this vsite and its contructing
1536 * atoms are all within the same cg, so no pbc is required.
1538 vsite_pbc_f[vsi] = -2;
1539 /* Check if constructing atoms are outside the vsite's cg */
1541 if (ftype == F_VSITEN) {
1542 nc3 = 3*iparams[ia[i]].vsiten.n;
1543 for(j=0; j<nc3; j+=3) {
1544 if (a2cg[ia[i+j+2]] != cg_v)
1545 vsite_pbc_f[vsi] = -1;
1548 for(a=1; a<nral; a++) {
1549 if (a2cg[ia[i+1+a]] != cg_v)
1550 vsite_pbc_f[vsi] = -1;
1553 if (vsite_pbc_f[vsi] == -1) {
1554 /* Check if this is the first processed atom of a vsite only cg */
1555 bViteOnlyCG_and_FirstAtom = TRUE;
1556 for(a=cgs->index[cg_v]; a<cgs->index[cg_v+1]; a++) {
1557 /* Non-vsites already have pbc set, so simply check for pbc_set */
1559 bViteOnlyCG_and_FirstAtom = FALSE;
1563 if (bViteOnlyCG_and_FirstAtom) {
1564 /* First processed atom of a vsite only charge group.
1565 * The pbc of the input coordinates to construct_vsites
1566 * should be preserved.
1568 vsite_pbc_f[vsi] = vsite;
1569 } else if (cg_v != a2cg[ia[1+i+1]]) {
1570 /* This vsite has a different charge group index
1571 * than it's first constructing atom
1572 * and the charge group has more than one atom,
1573 * search for the first normal particle
1574 * or vsite that already had its pbc defined.
1575 * If nothing is found, use full pbc for this vsite.
1577 for(a=cgs->index[cg_v]; a<cgs->index[cg_v+1]; a++) {
1578 if (a != vsite && pbc_set[a]) {
1579 vsite_pbc_f[vsi] = a;
1581 fprintf(debug,"vsite %d match pbc with atom %d\n",
1587 fprintf(debug,"vsite atom %d cg %d - %d pbc atom %d\n",
1588 vsite+1,cgs->index[cg_v]+1,cgs->index[cg_v+1],
1589 vsite_pbc_f[vsi]+1);
1592 if (ftype == F_VSITEN) {
1593 /* The other entries in vsite_pbc_f are not used for center vsites */
1599 /* This vsite now has its pbc defined */
1610 gmx_vsite_t *init_vsite(gmx_mtop_t *mtop,t_commrec *cr)
1616 gmx_moltype_t *molt;
1618 /* check if there are vsites */
1620 for(i=0; i<F_NRE; i++) {
1621 if (interaction_function[i].flags & IF_VSITE) {
1622 nvsite += gmx_mtop_ftype_count(mtop,i);
1632 vsite->n_intercg_vsite = count_intercg_vsite(mtop);
1634 if (vsite->n_intercg_vsite > 0 && DOMAINDECOMP(cr)) {
1635 vsite->nvsite_pbc_molt = mtop->nmoltype;
1636 snew(vsite->vsite_pbc_molt,vsite->nvsite_pbc_molt);
1637 for(mt=0; mt<mtop->nmoltype; mt++) {
1638 molt = &mtop->moltype[mt];
1639 /* Make an atom to charge group index */
1640 a2cg = atom2cg(&molt->cgs);
1641 vsite->vsite_pbc_molt[mt] = get_vsite_pbc(mtop->ffparams.iparams,
1643 molt->atoms.atom,NULL,
1648 snew(vsite->vsite_pbc_loc_nalloc,F_VSITEN-F_VSITE2+1);
1649 snew(vsite->vsite_pbc_loc ,F_VSITEN-F_VSITE2+1);
1656 void set_vsite_top(gmx_vsite_t *vsite,gmx_localtop_t *top,t_mdatoms *md,
1661 /* Make an atom to charge group index */
1662 a2cg = atom2cg(&top->cgs);
1664 if (vsite->n_intercg_vsite > 0) {
1665 vsite->vsite_pbc_loc = get_vsite_pbc(top->idef.iparams,
1666 top->idef.il,NULL,md,
1669 if (PARTDECOMP(cr)) {
1670 snew(vsite->vsitecomm,1);
1671 vsite->bPDvsitecomm =
1672 setup_parallel_vsites(&(top->idef),cr,vsite->vsitecomm);