?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
A freeze group can now be allowed to move rigidly in some dimension by using "freezed...
[gromacs/rigid-bodies.git] / src / mdlib / vsite.c
blobb2ce30d94d868429e9e307222feb7aa9bc1d651a
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * GROwing Monsters And Cloning Shrimps
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <stdio.h>
40 #include "typedefs.h"
41 #include "vsite.h"
42 #include "macros.h"
43 #include "smalloc.h"
44 #include "nrnb.h"
45 #include "vec.h"
46 #include "mvdata.h"
47 #include "network.h"
48 #include "mshift.h"
49 #include "pbc.h"
50 #include "domdec.h"
51 #include "partdec.h"
52 #include "mtop_util.h"
53
54 /* Routines to send/recieve coordinates and force
55 * of constructing atoms.
56 */
57
58 static void move_construct_x(t_comm_vsites *vsitecomm, rvec x[], t_commrec *cr)
59 {
60 rvec *sendbuf;
61 rvec *recvbuf;
62 int i,ia;
63
64 sendbuf = vsitecomm->send_buf;
65 recvbuf = vsitecomm->recv_buf;
66
67
68 /* Prepare pulse left by copying to send buffer */
69 for(i=0;i<vsitecomm->left_export_nconstruct;i++)
70 {
71 ia = vsitecomm->left_export_construct[i];
72 copy_rvec(x[ia],sendbuf[i]);
73 }
74
75 /* Pulse coordinates left */
76 gmx_tx_rx_real(cr,GMX_LEFT,(real *)sendbuf,3*vsitecomm->left_export_nconstruct,GMX_RIGHT,(real *)recvbuf,3*vsitecomm->right_import_nconstruct);
77
78 /* Copy from receive buffer to coordinate array */
79 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
80 {
81 ia = vsitecomm->right_import_construct[i];
82 copy_rvec(recvbuf[i],x[ia]);
83 }
84
85 /* Prepare pulse right by copying to send buffer */
86 for(i=0;i<vsitecomm->right_export_nconstruct;i++)
87 {
88 ia = vsitecomm->right_export_construct[i];
89 copy_rvec(x[ia],sendbuf[i]);
90 }
91
92 /* Pulse coordinates right */
93 gmx_tx_rx_real(cr,GMX_RIGHT,(real *)sendbuf,3*vsitecomm->right_export_nconstruct,GMX_LEFT,(real *)recvbuf,3*vsitecomm->left_import_nconstruct);
94
95 /* Copy from receive buffer to coordinate array */
96 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
97 {
98 ia = vsitecomm->left_import_construct[i];
99 copy_rvec(recvbuf[i],x[ia]);
100 }
101 }
102
103
104 static void move_construct_f(t_comm_vsites *vsitecomm, rvec f[], t_commrec *cr)
105 {
106 rvec *sendbuf;
107 rvec *recvbuf;
108 int i,ia;
109
110 sendbuf = vsitecomm->send_buf;
111 recvbuf = vsitecomm->recv_buf;
112
113 /* Prepare pulse right by copying to send buffer */
114 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
115 {
116 ia = vsitecomm->right_import_construct[i];
117 copy_rvec(f[ia],sendbuf[i]);
118 clear_rvec(f[ia]); /* Zero it here after moving, just to simplify debug book-keeping... */
119 }
120
121 /* Pulse forces right */
122 gmx_tx_rx_real(cr,GMX_RIGHT,(real *)sendbuf,3*vsitecomm->right_import_nconstruct,GMX_LEFT,(real *)recvbuf,3*vsitecomm->left_export_nconstruct);
123
124 /* Copy from receive buffer to coordinate array */
125 for(i=0;i<vsitecomm->left_export_nconstruct;i++)
126 {
127 ia = vsitecomm->left_export_construct[i];
128 rvec_inc(f[ia],recvbuf[i]);
129 }
130
131 /* Prepare pulse left by copying to send buffer */
132 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
133 {
134 ia = vsitecomm->left_import_construct[i];
135 copy_rvec(f[ia],sendbuf[i]);
136 clear_rvec(f[ia]); /* Zero it here after moving, just to simplify debug book-keeping... */
137 }
138
139 /* Pulse coordinates left */
140 gmx_tx_rx_real(cr,GMX_LEFT,(real *)sendbuf,3*vsitecomm->left_import_nconstruct,GMX_RIGHT,(real *)recvbuf,3*vsitecomm->right_export_nconstruct);
141
142 /* Copy from receive buffer to coordinate array */
143 for(i=0;i<vsitecomm->right_export_nconstruct;i++)
144 {
145 ia = vsitecomm->right_export_construct[i];
146 rvec_inc(f[ia],recvbuf[i]);
147 }
148
149 /* All forces are now on the home processors */
150 }
151
152
153 static void
154 pd_clear_nonlocal_constructs(t_comm_vsites *vsitecomm, rvec f[])
155 {
156 int i,ia;
157
158 for(i=0;i<vsitecomm->left_import_nconstruct;i++)
159 {
160 ia = vsitecomm->left_import_construct[i];
161 clear_rvec(f[ia]);
162 }
163 for(i=0;i<vsitecomm->right_import_nconstruct;i++)
164 {
165 ia = vsitecomm->right_import_construct[i];
166 clear_rvec(f[ia]);
167 }
168 }
169
170
171
172 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
173 {
174 if (pbc) {
175 return pbc_dx_aiuc(pbc,xi,xj,dx);
176 }
177 else {
178 rvec_sub(xi,xj,dx);
179 return CENTRAL;
180 }
181 }
182
183 /* Vsite construction routines */
184
185 static void constr_vsite2(rvec xi,rvec xj,rvec x,real a,t_pbc *pbc)
186 {
187 real b;
188 rvec dx;
189
190 b=1.0-a;
191 /* 1 flop */
192
193 if (pbc) {
194 pbc_dx_aiuc(pbc,xj,xi,dx);
195 x[XX] = xi[XX] + a*dx[XX];
196 x[YY] = xi[YY] + a*dx[YY];
197 x[ZZ] = xi[ZZ] + a*dx[ZZ];
198 } else {
199 x[XX] = b*xi[XX] + a*xj[XX];
200 x[YY] = b*xi[YY] + a*xj[YY];
201 x[ZZ] = b*xi[ZZ] + a*xj[ZZ];
202 /* 9 Flops */
203 }
204
205 /* TOTAL: 10 flops */
206 }
207
208 static void constr_vsite3(rvec xi,rvec xj,rvec xk,rvec x,real a,real b,
209 t_pbc *pbc)
210 {
211 real c;
212 rvec dxj,dxk;
213
214 c=1.0-a-b;
215 /* 2 flops */
216
217 if (pbc) {
218 pbc_dx_aiuc(pbc,xj,xi,dxj);
219 pbc_dx_aiuc(pbc,xk,xi,dxk);
220 x[XX] = xi[XX] + a*dxj[XX] + b*dxk[XX];
221 x[YY] = xi[YY] + a*dxj[YY] + b*dxk[YY];
222 x[ZZ] = xi[ZZ] + a*dxj[ZZ] + b*dxk[ZZ];
223 } else {
224 x[XX] = c*xi[XX] + a*xj[XX] + b*xk[XX];
225 x[YY] = c*xi[YY] + a*xj[YY] + b*xk[YY];
226 x[ZZ] = c*xi[ZZ] + a*xj[ZZ] + b*xk[ZZ];
227 /* 15 Flops */
228 }
229
230 /* TOTAL: 17 flops */
231 }
232
233 static void constr_vsite3FD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b,
234 t_pbc *pbc)
235 {
236 rvec xij,xjk,temp;
237 real c;
238
239 pbc_rvec_sub(pbc,xj,xi,xij);
240 pbc_rvec_sub(pbc,xk,xj,xjk);
241 /* 6 flops */
242
243 /* temp goes from i to a point on the line jk */
244 temp[XX] = xij[XX] + a*xjk[XX];
245 temp[YY] = xij[YY] + a*xjk[YY];
246 temp[ZZ] = xij[ZZ] + a*xjk[ZZ];
247 /* 6 flops */
248
249 c=b*gmx_invsqrt(iprod(temp,temp));
250 /* 6 + 10 flops */
251
252 x[XX] = xi[XX] + c*temp[XX];
253 x[YY] = xi[YY] + c*temp[YY];
254 x[ZZ] = xi[ZZ] + c*temp[ZZ];
255 /* 6 Flops */
256
257 /* TOTAL: 34 flops */
258 }
259
260 static void constr_vsite3FAD(rvec xi,rvec xj,rvec xk,rvec x,real a,real b, t_pbc *pbc)
261 {
262 rvec xij,xjk,xp;
263 real a1,b1,c1,invdij;
264
265 pbc_rvec_sub(pbc,xj,xi,xij);
266 pbc_rvec_sub(pbc,xk,xj,xjk);
267 /* 6 flops */
268
269 invdij = gmx_invsqrt(iprod(xij,xij));
270 c1 = invdij * invdij * iprod(xij,xjk);
271 xp[XX] = xjk[XX] - c1*xij[XX];
272 xp[YY] = xjk[YY] - c1*xij[YY];
273 xp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
274 a1 = a*invdij;
275 b1 = b*gmx_invsqrt(iprod(xp,xp));
276 /* 45 */
277
278 x[XX] = xi[XX] + a1*xij[XX] + b1*xp[XX];
279 x[YY] = xi[YY] + a1*xij[YY] + b1*xp[YY];
280 x[ZZ] = xi[ZZ] + a1*xij[ZZ] + b1*xp[ZZ];
281 /* 12 Flops */
282
283 /* TOTAL: 63 flops */
284 }
285
286 static void constr_vsite3OUT(rvec xi,rvec xj,rvec xk,rvec x,
287 real a,real b,real c,t_pbc *pbc)
288 {
289 rvec xij,xik,temp;
290
291 pbc_rvec_sub(pbc,xj,xi,xij);
292 pbc_rvec_sub(pbc,xk,xi,xik);
293 cprod(xij,xik,temp);
294 /* 15 Flops */
295
296 x[XX] = xi[XX] + a*xij[XX] + b*xik[XX] + c*temp[XX];
297 x[YY] = xi[YY] + a*xij[YY] + b*xik[YY] + c*temp[YY];
298 x[ZZ] = xi[ZZ] + a*xij[ZZ] + b*xik[ZZ] + c*temp[ZZ];
299 /* 18 Flops */
300
301 /* TOTAL: 33 flops */
302 }
303
304 static void constr_vsite4FD(rvec xi,rvec xj,rvec xk,rvec xl,rvec x,
305 real a,real b,real c,t_pbc *pbc)
306 {
307 rvec xij,xjk,xjl,temp;
308 real d;
309
310 pbc_rvec_sub(pbc,xj,xi,xij);
311 pbc_rvec_sub(pbc,xk,xj,xjk);
312 pbc_rvec_sub(pbc,xl,xj,xjl);
313 /* 9 flops */
314
315 /* temp goes from i to a point on the plane jkl */
316 temp[XX] = xij[XX] + a*xjk[XX] + b*xjl[XX];
317 temp[YY] = xij[YY] + a*xjk[YY] + b*xjl[YY];
318 temp[ZZ] = xij[ZZ] + a*xjk[ZZ] + b*xjl[ZZ];
319 /* 12 flops */
320
321 d=c*gmx_invsqrt(iprod(temp,temp));
322 /* 6 + 10 flops */
323
324 x[XX] = xi[XX] + d*temp[XX];
325 x[YY] = xi[YY] + d*temp[YY];
326 x[ZZ] = xi[ZZ] + d*temp[ZZ];
327 /* 6 Flops */
328
329 /* TOTAL: 43 flops */
330 }
331
332
333 static void constr_vsite4FDN(rvec xi,rvec xj,rvec xk,rvec xl,rvec x,
334 real a,real b,real c,t_pbc *pbc)
335 {
336 rvec xij,xik,xil,ra,rb,rja,rjb,rm;
337 real d;
338
339 pbc_rvec_sub(pbc,xj,xi,xij);
340 pbc_rvec_sub(pbc,xk,xi,xik);
341 pbc_rvec_sub(pbc,xl,xi,xil);
342 /* 9 flops */
343
344 ra[XX] = a*xik[XX];
345 ra[YY] = a*xik[YY];
346 ra[ZZ] = a*xik[ZZ];
347
348 rb[XX] = b*xil[XX];
349 rb[YY] = b*xil[YY];
350 rb[ZZ] = b*xil[ZZ];
351
352 /* 6 flops */
353
354 rvec_sub(ra,xij,rja);
355 rvec_sub(rb,xij,rjb);
356 /* 6 flops */
357
358 cprod(rja,rjb,rm);
359 /* 9 flops */
360
361 d=c*gmx_invsqrt(norm2(rm));
362 /* 5+5+1 flops */
363
364 x[XX] = xi[XX] + d*rm[XX];
365 x[YY] = xi[YY] + d*rm[YY];
366 x[ZZ] = xi[ZZ] + d*rm[ZZ];
367 /* 6 Flops */
368
369 /* TOTAL: 47 flops */
370 }
371
372
373 static int constr_vsiten(t_iatom *ia, t_iparams ip[],
374 rvec *x, t_pbc *pbc)
375 {
376 rvec xs,x1,dx;
377 dvec dsum;
378 int n3,av,ai,i;
379 real a;
380
381 n3 = 3*ip[ia[0]].vsiten.n;
382 av = ia[1];
383 ai = ia[2];
384 copy_rvec(x[ai],x1);
385 clear_dvec(dsum);
386 for(i=3; i<n3; i+=3) {
387 ai = ia[i+2];
388 a = ip[ia[i]].vsiten.a;
389 if (pbc) {
390 pbc_dx_aiuc(pbc,x[ai],x1,dx);
391 } else {
392 rvec_sub(x[ai],x1,dx);
393 }
394 dsum[XX] += a*dx[XX];
395 dsum[YY] += a*dx[YY];
396 dsum[ZZ] += a*dx[ZZ];
397 /* 9 Flops */
398 }
399
400 x[av][XX] = x1[XX] + dsum[XX];
401 x[av][YY] = x1[YY] + dsum[YY];
402 x[av][ZZ] = x1[ZZ] + dsum[ZZ];
403
404 return n3;
405 }
406
407
408 void construct_vsites(FILE *log,gmx_vsite_t *vsite,
409 rvec x[],t_nrnb *nrnb,
410 real dt,rvec *v,
411 t_iparams ip[],t_ilist ilist[],
412 int ePBC,gmx_bool bMolPBC,t_graph *graph,
413 t_commrec *cr,matrix box)
414 {
415 rvec xpbc,xv,vv,dx;
416 real a1,b1,c1,inv_dt;
417 int i,inc,ii,nra,nr,tp,ftype;
418 t_iatom avsite,ai,aj,ak,al,pbc_atom;
419 t_iatom *ia;
420 t_pbc pbc,*pbc_null,*pbc_null2;
421 gmx_bool bDomDec;
422 int *vsite_pbc,ishift;
423 rvec reftmp,vtmp,rtmp;
424
425 bDomDec = cr && DOMAINDECOMP(cr);
426
427 /* We only need to do pbc when we have inter-cg vsites */
428 if (ePBC != epbcNONE && (bDomDec || bMolPBC) && vsite->n_intercg_vsite) {
429 /* This is wasting some CPU time as we now do this multiple times
430 * per MD step. But how often do we have vsites with full pbc?
431 */
432 pbc_null = set_pbc_dd(&pbc,ePBC,cr!=NULL ? cr->dd : NULL,FALSE,box);
433 } else {
434 pbc_null = NULL;
435 }
436
437 if (cr) {
438 if (bDomDec) {
439 dd_move_x_vsites(cr->dd,box,x);
440 } else if (vsite->bPDvsitecomm) {
441 /* I'm not sure whether the periodicity and shift are guaranteed
442 * to be consistent between different nodes when running e.g. polymers
443 * in parallel. In this special case we thus unshift/shift,
444 * but only when necessary. This is to make sure the coordinates
445 * we move don't end up a box away...
446 */
447 if (graph)
448 unshift_self(graph,box,x);
449
450 move_construct_x(vsite->vsitecomm,x,cr);
451
452 if (graph)
453 shift_self(graph,box,x);
454 }
455 }
456
457 if (v) {
458 inv_dt = 1.0/dt;
459 } else {
460 inv_dt = 1.0;
461 }
462
463 pbc_null2 = NULL;
464 for(ftype=0; (ftype<F_NRE); ftype++) {
465 if (interaction_function[ftype].flags & IF_VSITE) {
466 nra = interaction_function[ftype].nratoms;
467 nr = ilist[ftype].nr;
468 ia = ilist[ftype].iatoms;
469
470 if (pbc_null) {
471 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
472 } else {
473 vsite_pbc = NULL;
474 }
475
476 for(i=0; (i<nr); ) {
477 tp = ia[0];
478 /*
479 if (ftype != idef->functype[tp])
480 gmx_incons("Function types for vsites wrong");
481 */
482
483 /* The vsite and constructing atoms */
484 avsite = ia[1];
485 ai = ia[2];
486 aj = ia[3];
487
488 /* Constants for constructing vsites */
489 a1 = ip[tp].vsite.a;
490 /* Check what kind of pbc we need to use */
491 if (vsite_pbc) {
492 pbc_atom = vsite_pbc[i/(1+nra)];
493 if (pbc_atom > -2) {
494 if (pbc_atom >= 0) {
495 /* We need to copy the coordinates here,
496 * single for single atom cg's pbc_atom is the vsite itself.
497 */
498 copy_rvec(x[pbc_atom],xpbc);
499 }
500 pbc_null2 = pbc_null;
501 } else {
502 pbc_null2 = NULL;
503 }
504 } else {
505 pbc_atom = -2;
506 }
507 /* Copy the old position */
508 copy_rvec(x[avsite],xv);
509
510 /* Construct the vsite depending on type */
511 inc = nra+1;
512 switch (ftype) {
513 case F_VSITE2:
514 constr_vsite2(x[ai],x[aj],x[avsite],a1,pbc_null2);
515 break;
516 case F_VSITE3:
517 ak = ia[4];
518 b1 = ip[tp].vsite.b;
519 constr_vsite3(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
520 break;
521 case F_VSITE3FD:
522 ak = ia[4];
523 b1 = ip[tp].vsite.b;
524 constr_vsite3FD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
525 break;
526 case F_VSITE3FAD:
527 ak = ia[4];
528 b1 = ip[tp].vsite.b;
529 constr_vsite3FAD(x[ai],x[aj],x[ak],x[avsite],a1,b1,pbc_null2);
530 break;
531 case F_VSITE3OUT:
532 ak = ia[4];
533 b1 = ip[tp].vsite.b;
534 c1 = ip[tp].vsite.c;
535 constr_vsite3OUT(x[ai],x[aj],x[ak],x[avsite],a1,b1,c1,pbc_null2);
536 break;
537 case F_VSITE4FD:
538 ak = ia[4];
539 al = ia[5];
540 b1 = ip[tp].vsite.b;
541 c1 = ip[tp].vsite.c;
542 constr_vsite4FD(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1,
543 pbc_null2);
544 break;
545 case F_VSITE4FDN:
546 ak = ia[4];
547 al = ia[5];
548 b1 = ip[tp].vsite.b;
549 c1 = ip[tp].vsite.c;
550 constr_vsite4FDN(x[ai],x[aj],x[ak],x[al],x[avsite],a1,b1,c1,
551 pbc_null2);
552 break;
553 case F_VSITEN:
554 inc = constr_vsiten(ia,ip,x,pbc_null2);
555 break;
556 default:
557 gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
558 ftype,__FILE__,__LINE__);
559 }
560
561 if (pbc_atom >= 0) {
562 /* Match the pbc of this vsite to the rest of its charge group */
563 ishift = pbc_dx_aiuc(pbc_null,x[avsite],xpbc,dx);
564 if (ishift != CENTRAL)
565 rvec_add(xpbc,dx,x[avsite]);
566 }
567 if (v) {
568 /* Calculate velocity of vsite... */
569 rvec_sub(x[avsite],xv,vv);
570 svmul(inv_dt,vv,v[avsite]);
571 }
572
573 /* Increment loop variables */
574 i += inc;
575 ia += inc;
576 }
577 }
578 }
579 }
580
581 static void spread_vsite2(t_iatom ia[],real a,
582 rvec x[],rvec f[],rvec fshift[],
583 t_pbc *pbc,t_graph *g)
584 {
585 rvec fi,fj,dx;
586 t_iatom av,ai,aj;
587 ivec di;
588 real b;
589 int siv,sij;
590
591 av = ia[1];
592 ai = ia[2];
593 aj = ia[3];
594
595 svmul(1-a,f[av],fi);
596 svmul( a,f[av],fj);
597 /* 7 flop */
598
599 rvec_inc(f[ai],fi);
600 rvec_inc(f[aj],fj);
601 /* 6 Flops */
602
603 if (g) {
604 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,av),di);
605 siv = IVEC2IS(di);
606 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di);
607 sij = IVEC2IS(di);
608 } else if (pbc) {
609 siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx);
610 sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
611 } else {
612 siv = CENTRAL;
613 sij = CENTRAL;
614 }
615
616 if (fshift && (siv != CENTRAL || sij != CENTRAL)) {
617 rvec_inc(fshift[siv],f[av]);
618 rvec_dec(fshift[CENTRAL],fi);
619 rvec_dec(fshift[sij],fj);
620 }
621
622 /* TOTAL: 13 flops */
623 }
624
625 void construct_vsites_mtop(FILE *log,gmx_vsite_t *vsite,
626 gmx_mtop_t *mtop,rvec x[])
627 {
628 int as,mb,mol;
629 gmx_molblock_t *molb;
630 gmx_moltype_t *molt;
631
632 as = 0;
633 for(mb=0; mb<mtop->nmolblock; mb++) {
634 molb = &mtop->molblock[mb];
635 molt = &mtop->moltype[molb->type];
636 for(mol=0; mol<molb->nmol; mol++) {
637 construct_vsites(log,vsite,x+as,NULL,0.0,NULL,
638 mtop->ffparams.iparams,molt->ilist,
639 epbcNONE,TRUE,NULL,NULL,NULL);
640 as += molt->atoms.nr;
641 }
642 }
643 }
644
645 static void spread_vsite3(t_iatom ia[],real a,real b,
646 rvec x[],rvec f[],rvec fshift[],
647 t_pbc *pbc,t_graph *g)
648 {
649 rvec fi,fj,fk,dx;
650 atom_id av,ai,aj,ak;
651 ivec di;
652 int siv,sij,sik;
653
654 av = ia[1];
655 ai = ia[2];
656 aj = ia[3];
657 ak = ia[4];
658
659 svmul(1-a-b,f[av],fi);
660 svmul( a,f[av],fj);
661 svmul( b,f[av],fk);
662 /* 11 flops */
663
664 rvec_inc(f[ai],fi);
665 rvec_inc(f[aj],fj);
666 rvec_inc(f[ak],fk);
667 /* 9 Flops */
668
669 if (g) {
670 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ia[1]),di);
671 siv = IVEC2IS(di);
672 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),di);
673 sij = IVEC2IS(di);
674 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),di);
675 sik = IVEC2IS(di);
676 } else if (pbc) {
677 siv = pbc_dx_aiuc(pbc,x[ai],x[av],dx);
678 sij = pbc_dx_aiuc(pbc,x[ai],x[aj],dx);
679 sik = pbc_dx_aiuc(pbc,x[ai],x[ak],dx);
680 } else {
681 siv = CENTRAL;
682 sij = CENTRAL;
683 sik = CENTRAL;
684 }
685
686 if (fshift && (siv!=CENTRAL || sij!=CENTRAL || sik!=CENTRAL)) {
687 rvec_inc(fshift[siv],f[av]);
688 rvec_dec(fshift[CENTRAL],fi);
689 rvec_dec(fshift[sij],fj);
690 rvec_dec(fshift[sik],fk);
691 }
692
693 /* TOTAL: 20 flops */
694 }
695
696 static void spread_vsite3FD(t_iatom ia[],real a,real b,
697 rvec x[],rvec f[],rvec fshift[],
698 t_pbc *pbc,t_graph *g)
699 {
700 real fx,fy,fz,c,invl,fproj,a1;
701 rvec xvi,xij,xjk,xix,fv,temp;
702 t_iatom av,ai,aj,ak;
703 int svi,sji,skj,d;
704 ivec di;
705
706 av = ia[1];
707 ai = ia[2];
708 aj = ia[3];
709 ak = ia[4];
710 copy_rvec(f[av],fv);
711
712 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
713 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
714 /* 6 flops */
715
716 /* xix goes from i to point x on the line jk */
717 xix[XX]=xij[XX]+a*xjk[XX];
718 xix[YY]=xij[YY]+a*xjk[YY];
719 xix[ZZ]=xij[ZZ]+a*xjk[ZZ];
720 /* 6 flops */
721
722 invl=gmx_invsqrt(iprod(xix,xix));
723 c=b*invl;
724 /* 4 + ?10? flops */
725
726 fproj=iprod(xix,fv)*invl*invl; /* = (xix . f)/(xix . xix) */
727
728 temp[XX]=c*(fv[XX]-fproj*xix[XX]);
729 temp[YY]=c*(fv[YY]-fproj*xix[YY]);
730 temp[ZZ]=c*(fv[ZZ]-fproj*xix[ZZ]);
731 /* 16 */
732
733 /* c is already calculated in constr_vsite3FD
734 storing c somewhere will save 26 flops! */
735
736 a1=1-a;
737 f[ai][XX] += fv[XX] - temp[XX];
738 f[ai][YY] += fv[YY] - temp[YY];
739 f[ai][ZZ] += fv[ZZ] - temp[ZZ];
740 f[aj][XX] += a1*temp[XX];
741 f[aj][YY] += a1*temp[YY];
742 f[aj][ZZ] += a1*temp[ZZ];
743 f[ak][XX] += a*temp[XX];
744 f[ak][YY] += a*temp[YY];
745 f[ak][ZZ] += a*temp[ZZ];
746 /* 19 Flops */
747
748 if (g) {
749 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
750 svi = IVEC2IS(di);
751 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
752 sji = IVEC2IS(di);
753 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
754 skj = IVEC2IS(di);
755 } else if (pbc) {
756 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
757 } else {
758 svi = CENTRAL;
759 }
760
761 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL)) {
762 rvec_dec(fshift[svi],fv);
763 fshift[CENTRAL][XX] += fv[XX] - (1 + a)*temp[XX];
764 fshift[CENTRAL][YY] += fv[YY] - (1 + a)*temp[YY];
765 fshift[CENTRAL][ZZ] += fv[ZZ] - (1 + a)*temp[ZZ];
766 fshift[ sji][XX] += temp[XX];
767 fshift[ sji][YY] += temp[YY];
768 fshift[ sji][ZZ] += temp[ZZ];
769 fshift[ skj][XX] += a*temp[XX];
770 fshift[ skj][YY] += a*temp[YY];
771 fshift[ skj][ZZ] += a*temp[ZZ];
772 }
773
774 /* TOTAL: 61 flops */
775 }
776
777 static void spread_vsite3FAD(t_iatom ia[],real a,real b,
778 rvec x[],rvec f[],rvec fshift[],
779 t_pbc *pbc,t_graph *g)
780 {
781 rvec xvi,xij,xjk,xperp,Fpij,Fppp,fv,f1,f2,f3;
782 real a1,b1,c1,c2,invdij,invdij2,invdp,fproj;
783 t_iatom av,ai,aj,ak;
784 int svi,sji,skj,d;
785 ivec di;
786
787 av = ia[1];
788 ai = ia[2];
789 aj = ia[3];
790 ak = ia[4];
791 copy_rvec(f[ia[1]],fv);
792
793 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
794 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
795 /* 6 flops */
796
797 invdij = gmx_invsqrt(iprod(xij,xij));
798 invdij2 = invdij * invdij;
799 c1 = iprod(xij,xjk) * invdij2;
800 xperp[XX] = xjk[XX] - c1*xij[XX];
801 xperp[YY] = xjk[YY] - c1*xij[YY];
802 xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ];
803 /* xperp in plane ijk, perp. to ij */
804 invdp = gmx_invsqrt(iprod(xperp,xperp));
805 a1 = a*invdij;
806 b1 = b*invdp;
807 /* 45 flops */
808
809 /* a1, b1 and c1 are already calculated in constr_vsite3FAD
810 storing them somewhere will save 45 flops! */
811
812 fproj=iprod(xij ,fv)*invdij2;
813 svmul(fproj, xij, Fpij); /* proj. f on xij */
814 svmul(iprod(xperp,fv)*invdp*invdp,xperp,Fppp); /* proj. f on xperp */
815 svmul(b1*fproj, xperp,f3);
816 /* 23 flops */
817
818 rvec_sub(fv,Fpij,f1); /* f1 = f - Fpij */
819 rvec_sub(f1,Fppp,f2); /* f2 = f - Fpij - Fppp */
820 for (d=0; (d<DIM); d++) {
821 f1[d]*=a1;
822 f2[d]*=b1;
823 }
824 /* 12 flops */
825
826 c2=1+c1;
827 f[ai][XX] += fv[XX] - f1[XX] + c1*f2[XX] + f3[XX];
828 f[ai][YY] += fv[YY] - f1[YY] + c1*f2[YY] + f3[YY];
829 f[ai][ZZ] += fv[ZZ] - f1[ZZ] + c1*f2[ZZ] + f3[ZZ];
830 f[aj][XX] += f1[XX] - c2*f2[XX] - f3[XX];
831 f[aj][YY] += f1[YY] - c2*f2[YY] - f3[YY];
832 f[aj][ZZ] += f1[ZZ] - c2*f2[ZZ] - f3[ZZ];
833 f[ak][XX] += f2[XX];
834 f[ak][YY] += f2[YY];
835 f[ak][ZZ] += f2[ZZ];
836 /* 30 Flops */
837
838 if (g) {
839 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
840 svi = IVEC2IS(di);
841 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
842 sji = IVEC2IS(di);
843 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
844 skj = IVEC2IS(di);
845 } else if (pbc) {
846 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
847 } else {
848 svi = CENTRAL;
849 }
850
851 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL)) {
852 rvec_dec(fshift[svi],fv);
853 fshift[CENTRAL][XX] += fv[XX] - f1[XX] - (1-c1)*f2[XX] + f3[XX];
854 fshift[CENTRAL][YY] += fv[YY] - f1[YY] - (1-c1)*f2[YY] + f3[YY];
855 fshift[CENTRAL][ZZ] += fv[ZZ] - f1[ZZ] - (1-c1)*f2[ZZ] + f3[ZZ];
856 fshift[ sji][XX] += f1[XX] - c1 *f2[XX] - f3[XX];
857 fshift[ sji][YY] += f1[YY] - c1 *f2[YY] - f3[YY];
858 fshift[ sji][ZZ] += f1[ZZ] - c1 *f2[ZZ] - f3[ZZ];
859 fshift[ skj][XX] += f2[XX];
860 fshift[ skj][YY] += f2[YY];
861 fshift[ skj][ZZ] += f2[ZZ];
862 }
863
864 /* TOTAL: 113 flops */
865 }
866
867 static void spread_vsite3OUT(t_iatom ia[],real a,real b,real c,
868 rvec x[],rvec f[],rvec fshift[],
869 t_pbc *pbc,t_graph *g)
870 {
871 rvec xvi,xij,xik,fv,fj,fk;
872 real cfx,cfy,cfz;
873 atom_id av,ai,aj,ak;
874 ivec di;
875 int svi,sji,ski;
876
877 av = ia[1];
878 ai = ia[2];
879 aj = ia[3];
880 ak = ia[4];
881
882 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
883 ski = pbc_rvec_sub(pbc,x[ak],x[ai],xik);
884 /* 6 Flops */
885
886 copy_rvec(f[av],fv);
887
888 cfx = c*fv[XX];
889 cfy = c*fv[YY];
890 cfz = c*fv[ZZ];
891 /* 3 Flops */
892
893 fj[XX] = a*fv[XX] - xik[ZZ]*cfy + xik[YY]*cfz;
894 fj[YY] = xik[ZZ]*cfx + a*fv[YY] - xik[XX]*cfz;
895 fj[ZZ] = -xik[YY]*cfx + xik[XX]*cfy + a*fv[ZZ];
896
897 fk[XX] = b*fv[XX] + xij[ZZ]*cfy - xij[YY]*cfz;
898 fk[YY] = -xij[ZZ]*cfx + b*fv[YY] + xij[XX]*cfz;
899 fk[ZZ] = xij[YY]*cfx - xij[XX]*cfy + b*fv[ZZ];
900 /* 30 Flops */
901
902 f[ai][XX] += fv[XX] - fj[XX] - fk[XX];
903 f[ai][YY] += fv[YY] - fj[YY] - fk[YY];
904 f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
905 rvec_inc(f[aj],fj);
906 rvec_inc(f[ak],fk);
907 /* 15 Flops */
908
909 if (g) {
910 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
911 svi = IVEC2IS(di);
912 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
913 sji = IVEC2IS(di);
914 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,ai),di);
915 ski = IVEC2IS(di);
916 } else if (pbc) {
917 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
918 } else {
919 svi = CENTRAL;
920 }
921
922 if (fshift && (svi!=CENTRAL || sji!=CENTRAL || ski!=CENTRAL)) {
923 rvec_dec(fshift[svi],fv);
924 fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX];
925 fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY];
926 fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ];
927 rvec_inc(fshift[sji],fj);
928 rvec_inc(fshift[ski],fk);
929 }
930
931 /* TOTAL: 54 flops */
932 }
933
934 static void spread_vsite4FD(t_iatom ia[],real a,real b,real c,
935 rvec x[],rvec f[],rvec fshift[],
936 t_pbc *pbc,t_graph *g)
937 {
938 real d,invl,fproj,a1;
939 rvec xvi,xij,xjk,xjl,xix,fv,temp;
940 atom_id av,ai,aj,ak,al;
941 ivec di;
942 int svi,sji,skj,slj,m;
943
944 av = ia[1];
945 ai = ia[2];
946 aj = ia[3];
947 ak = ia[4];
948 al = ia[5];
949
950 sji = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
951 skj = pbc_rvec_sub(pbc,x[ak],x[aj],xjk);
952 slj = pbc_rvec_sub(pbc,x[al],x[aj],xjl);
953 /* 9 flops */
954
955 /* xix goes from i to point x on the plane jkl */
956 for(m=0; m<DIM; m++)
957 xix[m] = xij[m] + a*xjk[m] + b*xjl[m];
958 /* 12 flops */
959
960 invl=gmx_invsqrt(iprod(xix,xix));
961 d=c*invl;
962 /* 4 + ?10? flops */
963
964 copy_rvec(f[av],fv);
965
966 fproj=iprod(xix,fv)*invl*invl; /* = (xix . f)/(xix . xix) */
967
968 for(m=0; m<DIM; m++)
969 temp[m] = d*(fv[m] - fproj*xix[m]);
970 /* 16 */
971
972 /* c is already calculated in constr_vsite3FD
973 storing c somewhere will save 35 flops! */
974
975 a1 = 1 - a - b;
976 for(m=0; m<DIM; m++) {
977 f[ai][m] += fv[m] - temp[m];
978 f[aj][m] += a1*temp[m];
979 f[ak][m] += a*temp[m];
980 f[al][m] += b*temp[m];
981 }
982 /* 26 Flops */
983
984 if (g) {
985 ivec_sub(SHIFT_IVEC(g,ia[1]),SHIFT_IVEC(g,ai),di);
986 svi = IVEC2IS(di);
987 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
988 sji = IVEC2IS(di);
989 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,aj),di);
990 skj = IVEC2IS(di);
991 ivec_sub(SHIFT_IVEC(g,al),SHIFT_IVEC(g,aj),di);
992 slj = IVEC2IS(di);
993 } else if (pbc) {
994 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
995 } else {
996 svi = CENTRAL;
997 }
998
999 if (fshift &&
1000 (svi!=CENTRAL || sji!=CENTRAL || skj!=CENTRAL || slj!=CENTRAL)) {
1001 rvec_dec(fshift[svi],fv);
1002 for(m=0; m<DIM; m++) {
1003 fshift[CENTRAL][m] += fv[m] - (1 + a + b)*temp[m];
1004 fshift[ sji][m] += temp[m];
1005 fshift[ skj][m] += a*temp[m];
1006 fshift[ slj][m] += b*temp[m];
1007 }
1008 }
1009
1010 /* TOTAL: 77 flops */
1011 }
1012
1013
1014 static void spread_vsite4FDN(t_iatom ia[],real a,real b,real c,
1015 rvec x[],rvec f[],rvec fshift[],
1016 t_pbc *pbc,t_graph *g)
1017 {
1018 rvec xvi,xij,xik,xil,ra,rb,rja,rjb,rab,rm,rt;
1019 rvec fv,fj,fk,fl;
1020 real invrm,denom;
1021 real cfx,cfy,cfz;
1022 ivec di;
1023 int av,ai,aj,ak,al;
1024 int svi,sij,sik,sil;
1025
1026 /* DEBUG: check atom indices */
1027 av = ia[1];
1028 ai = ia[2];
1029 aj = ia[3];
1030 ak = ia[4];
1031 al = ia[5];
1032
1033 copy_rvec(f[av],fv);
1034
1035 sij = pbc_rvec_sub(pbc,x[aj],x[ai],xij);
1036 sik = pbc_rvec_sub(pbc,x[ak],x[ai],xik);
1037 sil = pbc_rvec_sub(pbc,x[al],x[ai],xil);
1038 /* 9 flops */
1039
1040 ra[XX] = a*xik[XX];
1041 ra[YY] = a*xik[YY];
1042 ra[ZZ] = a*xik[ZZ];
1043
1044 rb[XX] = b*xil[XX];
1045 rb[YY] = b*xil[YY];
1046 rb[ZZ] = b*xil[ZZ];
1047
1048 /* 6 flops */
1049
1050 rvec_sub(ra,xij,rja);
1051 rvec_sub(rb,xij,rjb);
1052 rvec_sub(rb,ra,rab);
1053 /* 9 flops */
1054
1055 cprod(rja,rjb,rm);
1056 /* 9 flops */
1057
1058 invrm=gmx_invsqrt(norm2(rm));
1059 denom=invrm*invrm;
1060 /* 5+5+2 flops */
1061
1062 cfx = c*invrm*fv[XX];
1063 cfy = c*invrm*fv[YY];
1064 cfz = c*invrm*fv[ZZ];
1065 /* 6 Flops */
1066
1067 cprod(rm,rab,rt);
1068 /* 9 flops */
1069
1070 rt[XX] *= denom;
1071 rt[YY] *= denom;
1072 rt[ZZ] *= denom;
1073 /* 3flops */
1074
1075 fj[XX] = ( -rm[XX]*rt[XX]) * cfx + ( rab[ZZ]-rm[YY]*rt[XX]) * cfy + (-rab[YY]-rm[ZZ]*rt[XX]) * cfz;
1076 fj[YY] = (-rab[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + ( rab[XX]-rm[ZZ]*rt[YY]) * cfz;
1077 fj[ZZ] = ( rab[YY]-rm[XX]*rt[ZZ]) * cfx + (-rab[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1078 /* 30 flops */
1079
1080 cprod(rjb,rm,rt);
1081 /* 9 flops */
1082
1083 rt[XX] *= denom*a;
1084 rt[YY] *= denom*a;
1085 rt[ZZ] *= denom*a;
1086 /* 3flops */
1087
1088 fk[XX] = ( -rm[XX]*rt[XX]) * cfx + (-a*rjb[ZZ]-rm[YY]*rt[XX]) * cfy + ( a*rjb[YY]-rm[ZZ]*rt[XX]) * cfz;
1089 fk[YY] = ( a*rjb[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + (-a*rjb[XX]-rm[ZZ]*rt[YY]) * cfz;
1090 fk[ZZ] = (-a*rjb[YY]-rm[XX]*rt[ZZ]) * cfx + ( a*rjb[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1091 /* 36 flops */
1092
1093 cprod(rm,rja,rt);
1094 /* 9 flops */
1095
1096 rt[XX] *= denom*b;
1097 rt[YY] *= denom*b;
1098 rt[ZZ] *= denom*b;
1099 /* 3flops */
1100
1101 fl[XX] = ( -rm[XX]*rt[XX]) * cfx + ( b*rja[ZZ]-rm[YY]*rt[XX]) * cfy + (-b*rja[YY]-rm[ZZ]*rt[XX]) * cfz;
1102 fl[YY] = (-b*rja[ZZ]-rm[XX]*rt[YY]) * cfx + ( -rm[YY]*rt[YY]) * cfy + ( b*rja[XX]-rm[ZZ]*rt[YY]) * cfz;
1103 fl[ZZ] = ( b*rja[YY]-rm[XX]*rt[ZZ]) * cfx + (-b*rja[XX]-rm[YY]*rt[ZZ]) * cfy + ( -rm[ZZ]*rt[ZZ]) * cfz;
1104 /* 36 flops */
1105
1106 f[ai][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1107 f[ai][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1108 f[ai][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1109 rvec_inc(f[aj],fj);
1110 rvec_inc(f[ak],fk);
1111 rvec_inc(f[al],fl);
1112 /* 21 flops */
1113
1114 if (g) {
1115 ivec_sub(SHIFT_IVEC(g,av),SHIFT_IVEC(g,ai),di);
1116 svi = IVEC2IS(di);
1117 ivec_sub(SHIFT_IVEC(g,aj),SHIFT_IVEC(g,ai),di);
1118 sij = IVEC2IS(di);
1119 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,ai),di);
1120 sik = IVEC2IS(di);
1121 ivec_sub(SHIFT_IVEC(g,al),SHIFT_IVEC(g,ai),di);
1122 sil = IVEC2IS(di);
1123 } else if (pbc) {
1124 svi = pbc_rvec_sub(pbc,x[av],x[ai],xvi);
1125 } else {
1126 svi = CENTRAL;
1127 }
1128
1129 if (fshift && (svi!=CENTRAL || sij!=CENTRAL || sik!=CENTRAL || sil!=CENTRAL)) {
1130 rvec_dec(fshift[svi],fv);
1131 fshift[CENTRAL][XX] += fv[XX] - fj[XX] - fk[XX] - fl[XX];
1132 fshift[CENTRAL][YY] += fv[YY] - fj[YY] - fk[YY] - fl[YY];
1133 fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ] - fk[ZZ] - fl[ZZ];
1134 rvec_inc(fshift[sij],fj);
1135 rvec_inc(fshift[sik],fk);
1136 rvec_inc(fshift[sil],fl);
1137 }
1138
1139 /* Total: 207 flops (Yuck!) */
1140 }
1141
1142
1143 static int spread_vsiten(t_iatom ia[],t_iparams ip[],
1144 rvec x[],rvec f[],rvec fshift[],
1145 t_pbc *pbc,t_graph *g)
1146 {
1147 rvec xv,dx,fi;
1148 int n3,av,i,ai;
1149 real a;
1150 ivec di;
1151 int siv;
1152
1153 n3 = 3*ip[ia[0]].vsiten.n;
1154 av = ia[1];
1155 copy_rvec(x[av],xv);
1156
1157 for(i=0; i<n3; i+=3) {
1158 ai = ia[i+2];
1159 if (g) {
1160 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,av),di);
1161 siv = IVEC2IS(di);
1162 } else if (pbc) {
1163 siv = pbc_dx_aiuc(pbc,x[ai],xv,dx);
1164 } else {
1165 siv = CENTRAL;
1166 }
1167 a = ip[ia[i]].vsiten.a;
1168 svmul(a,f[av],fi);
1169 rvec_inc(f[ai],fi);
1170 if (fshift && siv != CENTRAL) {
1171 rvec_inc(fshift[siv],fi);
1172 rvec_dec(fshift[CENTRAL],fi);
1173 }
1174 /* 6 Flops */
1175 }
1176
1177 return n3;
1178 }
1179
1180
1181 void spread_vsite_f(FILE *log,gmx_vsite_t *vsite,
1182 rvec x[],rvec f[],rvec *fshift,
1183 t_nrnb *nrnb,t_idef *idef,
1184 int ePBC,gmx_bool bMolPBC,t_graph *g,matrix box,
1185 t_commrec *cr)
1186 {
1187 real a1,b1,c1;
1188 int i,inc,m,nra,nr,tp,ftype;
1189 int nd2,nd3,nd3FD,nd3FAD,nd3OUT,nd4FD,nd4FDN,ndN;
1190 t_iatom *ia;
1191 t_iparams *ip;
1192 t_pbc pbc,*pbc_null,*pbc_null2;
1193 int *vsite_pbc;
1194
1195 /* We only need to do pbc when we have inter-cg vsites */
1196 if ((DOMAINDECOMP(cr) || bMolPBC) && vsite->n_intercg_vsite) {
1197 /* This is wasting some CPU time as we now do this multiple times
1198 * per MD step. But how often do we have vsites with full pbc?
1199 */
1200 pbc_null = set_pbc_dd(&pbc,ePBC,cr->dd,FALSE,box);
1201 } else {
1202 pbc_null = NULL;
1203 }
1204
1205 if (DOMAINDECOMP(cr))
1206 {
1207 dd_clear_f_vsites(cr->dd,f);
1208 }
1209 else if (PARTDECOMP(cr) && vsite->vsitecomm != NULL)
1210 {
1211 pd_clear_nonlocal_constructs(vsite->vsitecomm,f);
1212 }
1213
1214 ip = idef->iparams;
1215
1216 nd2 = 0;
1217 nd3 = 0;
1218 nd3FD = 0;
1219 nd3FAD = 0;
1220 nd3OUT = 0;
1221 nd4FD = 0;
1222 nd4FDN = 0;
1223 ndN = 0;
1224
1225 /* this loop goes backwards to be able to build *
1226 * higher type vsites from lower types */
1227 pbc_null2 = NULL;
1228 for(ftype=F_NRE-1; (ftype>=0); ftype--) {
1229 if (interaction_function[ftype].flags & IF_VSITE) {
1230 nra = interaction_function[ftype].nratoms;
1231 nr = idef->il[ftype].nr;
1232 ia = idef->il[ftype].iatoms;
1233
1234 if (pbc_null) {
1235 vsite_pbc = vsite->vsite_pbc_loc[ftype-F_VSITE2];
1236 } else {
1237 vsite_pbc = NULL;
1238 }
1239
1240 for(i=0; (i<nr); ) {
1241 /* Check if we need to apply pbc for this vsite */
1242 if (vsite_pbc) {
1243 if (vsite_pbc[i/(1+nra)] > -2)
1244 pbc_null2 = pbc_null;
1245 else
1246 pbc_null2 = NULL;
1247 }
1248
1249 tp = ia[0];
1250 if (ftype != idef->functype[tp])
1251 gmx_incons("Functiontypes for vsites wrong");
1252
1253 /* Constants for constructing */
1254 a1 = ip[tp].vsite.a;
1255 /* Construct the vsite depending on type */
1256 inc = nra+1;
1257 switch (ftype) {
1258 case F_VSITE2:
1259 spread_vsite2(ia,a1,x,f,fshift,pbc_null2,g);
1260 nd2++;
1261 break;
1262 case F_VSITE3:
1263 b1 = ip[tp].vsite.b;
1264 spread_vsite3(ia,a1,b1,x,f,fshift,pbc_null2,g);
1265 nd3++;
1266 break;
1267 case F_VSITE3FD:
1268 b1 = ip[tp].vsite.b;
1269 spread_vsite3FD(ia,a1,b1,x,f,fshift,pbc_null2,g);
1270 nd3FD++;
1271 break;
1272 case F_VSITE3FAD:
1273 b1 = ip[tp].vsite.b;
1274 spread_vsite3FAD(ia,a1,b1,x,f,fshift,pbc_null2,g);
1275 nd3FAD++;
1276 break;
1277 case F_VSITE3OUT:
1278 b1 = ip[tp].vsite.b;
1279 c1 = ip[tp].vsite.c;
1280 spread_vsite3OUT(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
1281 nd3OUT++;
1282 break;
1283 case F_VSITE4FD:
1284 b1 = ip[tp].vsite.b;
1285 c1 = ip[tp].vsite.c;
1286 spread_vsite4FD(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
1287 nd4FD++;
1288 break;
1289 case F_VSITE4FDN:
1290 b1 = ip[tp].vsite.b;
1291 c1 = ip[tp].vsite.c;
1292 spread_vsite4FDN(ia,a1,b1,c1,x,f,fshift,pbc_null2,g);
1293 nd4FDN++;
1294 break;
1295 case F_VSITEN:
1296 inc = spread_vsiten(ia,ip,x,f,fshift,pbc_null2,g);
1297 ndN += inc;
1298 break;
1299 default:
1300 gmx_fatal(FARGS,"No such vsite type %d in %s, line %d",
1301 ftype,__FILE__,__LINE__);
1302 }
1303 clear_rvec(f[ia[1]]);
1304
1305 /* Increment loop variables */
1306 i += inc;
1307 ia += inc;
1308 }
1309 }
1310 }
1311
1312 inc_nrnb(nrnb,eNR_VSITE2, nd2 );
1313 inc_nrnb(nrnb,eNR_VSITE3, nd3 );
1314 inc_nrnb(nrnb,eNR_VSITE3FD, nd3FD );
1315 inc_nrnb(nrnb,eNR_VSITE3FAD,nd3FAD );
1316 inc_nrnb(nrnb,eNR_VSITE3OUT,nd3OUT );
1317 inc_nrnb(nrnb,eNR_VSITE4FD, nd4FD );
1318 inc_nrnb(nrnb,eNR_VSITE4FDN,nd4FDN );
1319 inc_nrnb(nrnb,eNR_VSITEN, ndN );
1320
1321 if (DOMAINDECOMP(cr)) {
1322 dd_move_f_vsites(cr->dd,f,fshift);
1323 } else if (vsite->bPDvsitecomm) {
1324 /* We only move forces here, and they are independent of shifts */
1325 move_construct_f(vsite->vsitecomm,f,cr);
1326 }
1327 }
1328
1329 static int *atom2cg(t_block *cgs)
1330 {
1331 int *a2cg,cg,i;
1332
1333 snew(a2cg,cgs->index[cgs->nr]);
1334 for(cg=0; cg<cgs->nr; cg++) {
1335 for(i=cgs->index[cg]; i<cgs->index[cg+1]; i++)
1336 a2cg[i] = cg;
1337 }
1338
1339 return a2cg;
1340 }
1341
1342 static int count_intercg_vsite(gmx_mtop_t *mtop)
1343 {
1344 int mb,mt,ftype,nral,i,cg,a;
1345 gmx_molblock_t *molb;
1346 gmx_moltype_t *molt;
1347 int *a2cg;
1348 t_ilist *il;
1349 t_iatom *ia;
1350 int n_intercg_vsite;
1351
1352 n_intercg_vsite = 0;
1353 for(mb=0; mb<mtop->nmolblock; mb++) {
1354 molb = &mtop->molblock[mb];
1355 molt = &mtop->moltype[molb->type];
1356 a2cg = atom2cg(&molt->cgs);
1357 for(ftype=0; ftype<F_NRE; ftype++) {
1358 if (interaction_function[ftype].flags & IF_VSITE) {
1359 nral = NRAL(ftype);
1360 il = &molt->ilist[ftype];
1361 ia = il->iatoms;
1362 for(i=0; i<il->nr; i+=1+nral) {
1363 cg = a2cg[ia[1+i]];
1364 for(a=1; a<nral; a++) {
1365 if (a2cg[ia[1+a]] != cg) {
1366 n_intercg_vsite += molb->nmol;
1367 break;
1368 }
1369 }
1370 }
1371 }
1372 }
1373 sfree(a2cg);
1374 }
1375
1376 return n_intercg_vsite;
1377 }
1378
1379 static int **get_vsite_pbc(t_iparams *iparams,t_ilist *ilist,
1380 t_atom *atom,t_mdatoms *md,
1381 t_block *cgs,int *a2cg)
1382 {
1383 int ftype,nral,i,j,vsi,vsite,cg_v,cg_c,a,nc3=0;
1384 t_ilist *il;
1385 t_iatom *ia;
1386 int **vsite_pbc,*vsite_pbc_f;
1387 char *pbc_set;
1388 gmx_bool bViteOnlyCG_and_FirstAtom;
1389
1390 /* Make an array that tells if the pbc of an atom is set */
1391 snew(pbc_set,cgs->index[cgs->nr]);
1392 /* PBC is set for all non vsites */
1393 for(a=0; a<cgs->index[cgs->nr]; a++) {
1394 if ((atom && atom[a].ptype != eptVSite) ||
1395 (md && md->ptype[a] != eptVSite)) {
1396 pbc_set[a] = 1;
1397 }
1398 }
1399
1400 snew(vsite_pbc,F_VSITEN-F_VSITE2+1);
1401
1402 for(ftype=0; ftype<F_NRE; ftype++) {
1403 if (interaction_function[ftype].flags & IF_VSITE) {
1404 nral = NRAL(ftype);
1405 il = &ilist[ftype];
1406 ia = il->iatoms;
1407
1408 snew(vsite_pbc[ftype-F_VSITE2],il->nr/(1+nral));
1409 vsite_pbc_f = vsite_pbc[ftype-F_VSITE2];
1410
1411 i = 0;
1412 while (i < il->nr) {
1413 vsi = i/(1+nral);
1414 vsite = ia[i+1];
1415 cg_v = a2cg[vsite];
1416 /* A value of -2 signals that this vsite and its contructing
1417 * atoms are all within the same cg, so no pbc is required.
1418 */
1419 vsite_pbc_f[vsi] = -2;
1420 /* Check if constructing atoms are outside the vsite's cg */
1421 nc3 = 0;
1422 if (ftype == F_VSITEN) {
1423 nc3 = 3*iparams[ia[i]].vsiten.n;
1424 for(j=0; j<nc3; j+=3) {
1425 if (a2cg[ia[i+j+2]] != cg_v)
1426 vsite_pbc_f[vsi] = -1;
1427 }
1428 } else {
1429 for(a=1; a<nral; a++) {
1430 if (a2cg[ia[i+1+a]] != cg_v)
1431 vsite_pbc_f[vsi] = -1;
1432 }
1433 }
1434 if (vsite_pbc_f[vsi] == -1) {
1435 /* Check if this is the first processed atom of a vsite only cg */
1436 bViteOnlyCG_and_FirstAtom = TRUE;
1437 for(a=cgs->index[cg_v]; a<cgs->index[cg_v+1]; a++) {
1438 /* Non-vsites already have pbc set, so simply check for pbc_set */
1439 if (pbc_set[a]) {
1440 bViteOnlyCG_and_FirstAtom = FALSE;
1441 break;
1442 }
1443 }
1444 if (bViteOnlyCG_and_FirstAtom) {
1445 /* First processed atom of a vsite only charge group.
1446 * The pbc of the input coordinates to construct_vsites
1447 * should be preserved.
1448 */
1449 vsite_pbc_f[vsi] = vsite;
1450 } else if (cg_v != a2cg[ia[1+i+1]]) {
1451 /* This vsite has a different charge group index
1452 * than it's first constructing atom
1453 * and the charge group has more than one atom,
1454 * search for the first normal particle
1455 * or vsite that already had its pbc defined.
1456 * If nothing is found, use full pbc for this vsite.
1457 */
1458 for(a=cgs->index[cg_v]; a<cgs->index[cg_v+1]; a++) {
1459 if (a != vsite && pbc_set[a]) {
1460 vsite_pbc_f[vsi] = a;
1461 if (gmx_debug_at)
1462 fprintf(debug,"vsite %d match pbc with atom %d\n",
1463 vsite+1,a+1);
1464 break;
1465 }
1466 }
1467 if (gmx_debug_at)
1468 fprintf(debug,"vsite atom %d cg %d - %d pbc atom %d\n",
1469 vsite+1,cgs->index[cg_v]+1,cgs->index[cg_v+1],
1470 vsite_pbc_f[vsi]+1);
1471 }
1472 }
1473 if (ftype == F_VSITEN) {
1474 /* The other entries in vsite_pbc_f are not used for center vsites */
1475 i += nc3;
1476 } else {
1477 i += 1+nral;
1478 }
1479
1480 /* This vsite now has its pbc defined */
1481 pbc_set[vsite] = 1;
1482 }
1483 }
1484 }
1485
1486 sfree(pbc_set);
1487
1488 return vsite_pbc;
1489 }
1490
1491 gmx_vsite_t *init_vsite(gmx_mtop_t *mtop,t_commrec *cr)
1492 {
1493 int nvsite,i;
1494 int *a2cg,cg;
1495 gmx_vsite_t *vsite;
1496 int mt;
1497 gmx_moltype_t *molt;
1498
1499 /* check if there are vsites */
1500 nvsite = 0;
1501 for(i=0; i<F_NRE; i++) {
1502 if (interaction_function[i].flags & IF_VSITE) {
1503 nvsite += gmx_mtop_ftype_count(mtop,i);
1504 }
1505 }
1506
1507 if (nvsite == 0) {
1508 return NULL;
1509 }
1510
1511 snew(vsite,1);
1512
1513 vsite->n_intercg_vsite = count_intercg_vsite(mtop);
1514
1515 if (vsite->n_intercg_vsite > 0 && DOMAINDECOMP(cr)) {
1516 vsite->nvsite_pbc_molt = mtop->nmoltype;
1517 snew(vsite->vsite_pbc_molt,vsite->nvsite_pbc_molt);
1518 for(mt=0; mt<mtop->nmoltype; mt++) {
1519 molt = &mtop->moltype[mt];
1520 /* Make an atom to charge group index */
1521 a2cg = atom2cg(&molt->cgs);
1522 vsite->vsite_pbc_molt[mt] = get_vsite_pbc(mtop->ffparams.iparams,
1523 molt->ilist,
1524 molt->atoms.atom,NULL,
1525 &molt->cgs,a2cg);
1526 sfree(a2cg);
1527 }
1528
1529 snew(vsite->vsite_pbc_loc_nalloc,F_VSITEN-F_VSITE2+1);
1530 snew(vsite->vsite_pbc_loc ,F_VSITEN-F_VSITE2+1);
1531 }
1532
1533
1534 return vsite;
1535 }
1536
1537 void set_vsite_top(gmx_vsite_t *vsite,gmx_localtop_t *top,t_mdatoms *md,
1538 t_commrec *cr)
1539 {
1540 int *a2cg;
1541
1542 /* Make an atom to charge group index */
1543 a2cg = atom2cg(&top->cgs);
1544
1545 if (vsite->n_intercg_vsite > 0) {
1546 vsite->vsite_pbc_loc = get_vsite_pbc(top->idef.iparams,
1547 top->idef.il,NULL,md,
1548 &top->cgs,a2cg);
1549
1550 if (PARTDECOMP(cr)) {
1551 snew(vsite->vsitecomm,1);
1552 vsite->bPDvsitecomm =
1553 setup_parallel_vsites(&(top->idef),cr,vsite->vsitecomm);
1554 }
1555 }
1556
1557 sfree(a2cg);
1558 }