?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Revised wording in pdb2gmx.c, hopefully clearer now.
[gromacs/rigid-bodies.git] / src / gmxlib / bondfree.c
blobf27df6d4790fb925c4366de6b3c8f1d67894a69f
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
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.
15
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.
20 *
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.
27 *
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.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROningen Mixture of Alchemy and Childrens' Stories
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include "physics.h"
42 #include "vec.h"
43 #include "maths.h"
44 #include "txtdump.h"
45 #include "bondf.h"
46 #include "smalloc.h"
47 #include "pbc.h"
48 #include "ns.h"
49 #include "macros.h"
50 #include "names.h"
51 #include "gmx_fatal.h"
52 #include "mshift.h"
53 #include "main.h"
54 #include "disre.h"
55 #include "orires.h"
56 #include "force.h"
57 #include "nonbonded.h"
58 #include "mdrun.h"
59
60 /* Find a better place for this? */
61 const int cmap_coeff_matrix[] = {
62 1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4 ,
63 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
64 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4 ,
65 0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
66 0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2 ,
67 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2 ,
68 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
69 0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
70 0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2 ,
71 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
72 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
73 0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2 ,
74 0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
75 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
76 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
77 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
78 };
79
80
81
82 int glatnr(int *global_atom_index,int i)
83 {
84 int atnr;
85
86 if (global_atom_index == NULL) {
87 atnr = i + 1;
88 } else {
89 atnr = global_atom_index[i] + 1;
90 }
91
92 return atnr;
93 }
94
95 static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
96 {
97 if (pbc) {
98 return pbc_dx_aiuc(pbc,xi,xj,dx);
99 }
100 else {
101 rvec_sub(xi,xj,dx);
102 return CENTRAL;
103 }
104 }
105
106 /*
107 * Morse potential bond by Frank Everdij
108 *
109 * Three parameters needed:
110 *
111 * b0 = equilibrium distance in nm
112 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
113 * cb = well depth in kJ/mol
114 *
115 * Note: the potential is referenced to be +cb at infinite separation
116 * and zero at the equilibrium distance!
117 */
118
119 real morse_bonds(int nbonds,
120 const t_iatom forceatoms[],const t_iparams forceparams[],
121 const rvec x[],rvec f[],rvec fshift[],
122 const t_pbc *pbc,const t_graph *g,
123 real lambda,real *dvdl,
124 const t_mdatoms *md,t_fcdata *fcd,
125 int *global_atom_index)
126 {
127 const real one=1.0;
128 const real two=2.0;
129 real dr,dr2,temp,omtemp,cbomtemp,fbond,vbond,fij,b0,be,cb,vtot;
130 rvec dx;
131 int i,m,ki,type,ai,aj;
132 ivec dt;
133
134 vtot = 0.0;
135 for(i=0; (i<nbonds); ) {
136 type = forceatoms[i++];
137 ai = forceatoms[i++];
138 aj = forceatoms[i++];
139
140 b0 = forceparams[type].morse.b0;
141 be = forceparams[type].morse.beta;
142 cb = forceparams[type].morse.cb;
143
144 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
145 dr2 = iprod(dx,dx); /* 5 */
146 dr = dr2*gmx_invsqrt(dr2); /* 10 */
147 temp = exp(-be*(dr-b0)); /* 12 */
148
149 if (temp == one)
150 continue;
151
152 omtemp = one-temp; /* 1 */
153 cbomtemp = cb*omtemp; /* 1 */
154 vbond = cbomtemp*omtemp; /* 1 */
155 fbond = -two*be*temp*cbomtemp*gmx_invsqrt(dr2); /* 9 */
156 vtot += vbond; /* 1 */
157
158 if (g) {
159 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
160 ki = IVEC2IS(dt);
161 }
162
163 for (m=0; (m<DIM); m++) { /* 15 */
164 fij=fbond*dx[m];
165 f[ai][m]+=fij;
166 f[aj][m]-=fij;
167 fshift[ki][m]+=fij;
168 fshift[CENTRAL][m]-=fij;
169 }
170 } /* 58 TOTAL */
171 return vtot;
172 }
173
174 real cubic_bonds(int nbonds,
175 const t_iatom forceatoms[],const t_iparams forceparams[],
176 const rvec x[],rvec f[],rvec fshift[],
177 const t_pbc *pbc,const t_graph *g,
178 real lambda,real *dvdl,
179 const t_mdatoms *md,t_fcdata *fcd,
180 int *global_atom_index)
181 {
182 const real three = 3.0;
183 const real two = 2.0;
184 real kb,b0,kcub;
185 real dr,dr2,dist,kdist,kdist2,fbond,vbond,fij,vtot;
186 rvec dx;
187 int i,m,ki,type,ai,aj;
188 ivec dt;
189
190 vtot = 0.0;
191 for(i=0; (i<nbonds); ) {
192 type = forceatoms[i++];
193 ai = forceatoms[i++];
194 aj = forceatoms[i++];
195
196 b0 = forceparams[type].cubic.b0;
197 kb = forceparams[type].cubic.kb;
198 kcub = forceparams[type].cubic.kcub;
199
200 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
201 dr2 = iprod(dx,dx); /* 5 */
202
203 if (dr2 == 0.0)
204 continue;
205
206 dr = dr2*gmx_invsqrt(dr2); /* 10 */
207 dist = dr-b0;
208 kdist = kb*dist;
209 kdist2 = kdist*dist;
210
211 vbond = kdist2 + kcub*kdist2*dist;
212 fbond = -(two*kdist + three*kdist2*kcub)/dr;
213
214 vtot += vbond; /* 21 */
215
216 if (g) {
217 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
218 ki=IVEC2IS(dt);
219 }
220 for (m=0; (m<DIM); m++) { /* 15 */
221 fij=fbond*dx[m];
222 f[ai][m]+=fij;
223 f[aj][m]-=fij;
224 fshift[ki][m]+=fij;
225 fshift[CENTRAL][m]-=fij;
226 }
227 } /* 54 TOTAL */
228 return vtot;
229 }
230
231 real FENE_bonds(int nbonds,
232 const t_iatom forceatoms[],const t_iparams forceparams[],
233 const rvec x[],rvec f[],rvec fshift[],
234 const t_pbc *pbc,const t_graph *g,
235 real lambda,real *dvdl,
236 const t_mdatoms *md,t_fcdata *fcd,
237 int *global_atom_index)
238 {
239 const real half=0.5;
240 const real one=1.0;
241 real bm,kb;
242 real dr,dr2,bm2,omdr2obm2,fbond,vbond,fij,vtot;
243 rvec dx;
244 int i,m,ki,type,ai,aj;
245 ivec dt;
246
247 vtot = 0.0;
248 for(i=0; (i<nbonds); ) {
249 type = forceatoms[i++];
250 ai = forceatoms[i++];
251 aj = forceatoms[i++];
252
253 bm = forceparams[type].fene.bm;
254 kb = forceparams[type].fene.kb;
255
256 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
257 dr2 = iprod(dx,dx); /* 5 */
258
259 if (dr2 == 0.0)
260 continue;
261
262 bm2 = bm*bm;
263
264 if (dr2 >= bm2)
265 gmx_fatal(FARGS,
266 "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d",
267 dr2,bm2,
268 glatnr(global_atom_index,ai),
269 glatnr(global_atom_index,aj));
270
271 omdr2obm2 = one - dr2/bm2;
272
273 vbond = -half*kb*bm2*log(omdr2obm2);
274 fbond = -kb/omdr2obm2;
275
276 vtot += vbond; /* 35 */
277
278 if (g) {
279 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
280 ki=IVEC2IS(dt);
281 }
282 for (m=0; (m<DIM); m++) { /* 15 */
283 fij=fbond*dx[m];
284 f[ai][m]+=fij;
285 f[aj][m]-=fij;
286 fshift[ki][m]+=fij;
287 fshift[CENTRAL][m]-=fij;
288 }
289 } /* 58 TOTAL */
290 return vtot;
291 }
292
293 real harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
294 real *V,real *F)
295 {
296 const real half=0.5;
297 real L1,kk,x0,dx,dx2;
298 real v,f,dvdl;
299
300 L1 = 1.0-lambda;
301 kk = L1*kA+lambda*kB;
302 x0 = L1*xA+lambda*xB;
303
304 dx = x-x0;
305 dx2 = dx*dx;
306
307 f = -kk*dx;
308 v = half*kk*dx2;
309 dvdl = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
310
311 *F = f;
312 *V = v;
313
314 return dvdl;
315
316 /* That was 19 flops */
317 }
318
319
320 real bonds(int nbonds,
321 const t_iatom forceatoms[],const t_iparams forceparams[],
322 const rvec x[],rvec f[],rvec fshift[],
323 const t_pbc *pbc,const t_graph *g,
324 real lambda,real *dvdlambda,
325 const t_mdatoms *md,t_fcdata *fcd,
326 int *global_atom_index)
327 {
328 int i,m,ki,ai,aj,type;
329 real dr,dr2,fbond,vbond,fij,vtot;
330 rvec dx;
331 ivec dt;
332
333 vtot = 0.0;
334 for(i=0; (i<nbonds); ) {
335 type = forceatoms[i++];
336 ai = forceatoms[i++];
337 aj = forceatoms[i++];
338
339 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
340 dr2 = iprod(dx,dx); /* 5 */
341 dr = dr2*gmx_invsqrt(dr2); /* 10 */
342
343 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
344 forceparams[type].harmonic.krB,
345 forceparams[type].harmonic.rA,
346 forceparams[type].harmonic.rB,
347 dr,lambda,&vbond,&fbond); /* 19 */
348
349 if (dr2 == 0.0)
350 continue;
351
352
353 vtot += vbond;/* 1*/
354 fbond *= gmx_invsqrt(dr2); /* 6 */
355 #ifdef DEBUG
356 if (debug)
357 fprintf(debug,"BONDS: dr = %10g vbond = %10g fbond = %10g\n",
358 dr,vbond,fbond);
359 #endif
360 if (g) {
361 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
362 ki=IVEC2IS(dt);
363 }
364 for (m=0; (m<DIM); m++) { /* 15 */
365 fij=fbond*dx[m];
366 f[ai][m]+=fij;
367 f[aj][m]-=fij;
368 fshift[ki][m]+=fij;
369 fshift[CENTRAL][m]-=fij;
370 }
371 } /* 59 TOTAL */
372 return vtot;
373 }
374
375 real restraint_bonds(int nbonds,
376 const t_iatom forceatoms[],const t_iparams forceparams[],
377 const rvec x[],rvec f[],rvec fshift[],
378 const t_pbc *pbc,const t_graph *g,
379 real lambda,real *dvdlambda,
380 const t_mdatoms *md,t_fcdata *fcd,
381 int *global_atom_index)
382 {
383 int i,m,ki,ai,aj,type;
384 real dr,dr2,fbond,vbond,fij,vtot;
385 real L1;
386 real low,dlow,up1,dup1,up2,dup2,k,dk;
387 real drh,drh2;
388 rvec dx;
389 ivec dt;
390
391 L1 = 1.0 - lambda;
392
393 vtot = 0.0;
394 for(i=0; (i<nbonds); )
395 {
396 type = forceatoms[i++];
397 ai = forceatoms[i++];
398 aj = forceatoms[i++];
399
400 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
401 dr2 = iprod(dx,dx); /* 5 */
402 dr = dr2*gmx_invsqrt(dr2); /* 10 */
403
404 low = L1*forceparams[type].restraint.lowA + lambda*forceparams[type].restraint.lowB;
405 dlow = -forceparams[type].restraint.lowA + forceparams[type].restraint.lowB;
406 up1 = L1*forceparams[type].restraint.up1A + lambda*forceparams[type].restraint.up1B;
407 dup1 = -forceparams[type].restraint.up1A + forceparams[type].restraint.up1B;
408 up2 = L1*forceparams[type].restraint.up2A + lambda*forceparams[type].restraint.up2B;
409 dup2 = -forceparams[type].restraint.up2A + forceparams[type].restraint.up2B;
410 k = L1*forceparams[type].restraint.kA + lambda*forceparams[type].restraint.kB;
411 dk = -forceparams[type].restraint.kA + forceparams[type].restraint.kB;
412 /* 24 */
413
414 if (dr < low)
415 {
416 drh = dr - low;
417 drh2 = drh*drh;
418 vbond = 0.5*k*drh2;
419 fbond = -k*drh;
420 *dvdlambda += 0.5*dk*drh2 - k*dlow*drh;
421 } /* 11 */
422 else if (dr <= up1)
423 {
424 vbond = 0;
425 fbond = 0;
426 }
427 else if (dr <= up2)
428 {
429 drh = dr - up1;
430 drh2 = drh*drh;
431 vbond = 0.5*k*drh2;
432 fbond = -k*drh;
433 *dvdlambda += 0.5*dk*drh2 - k*dup1*drh;
434 } /* 11 */
435 else
436 {
437 drh = dr - up2;
438 vbond = k*(up2 - up1)*(0.5*(up2 - up1) + drh);
439 fbond = -k*(up2 - up1);
440 *dvdlambda += dk*(up2 - up1)*(0.5*(up2 - up1) + drh)
441 + k*(dup2 - dup1)*(up2 - up1 + drh)
442 - k*(up2 - up1)*dup2;
443 }
444
445 if (dr2 == 0.0)
446 continue;
447
448 vtot += vbond;/* 1*/
449 fbond *= gmx_invsqrt(dr2); /* 6 */
450 #ifdef DEBUG
451 if (debug)
452 fprintf(debug,"BONDS: dr = %10g vbond = %10g fbond = %10g\n",
453 dr,vbond,fbond);
454 #endif
455 if (g) {
456 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
457 ki=IVEC2IS(dt);
458 }
459 for (m=0; (m<DIM); m++) { /* 15 */
460 fij=fbond*dx[m];
461 f[ai][m]+=fij;
462 f[aj][m]-=fij;
463 fshift[ki][m]+=fij;
464 fshift[CENTRAL][m]-=fij;
465 }
466 } /* 59 TOTAL */
467
468 return vtot;
469 }
470
471 real polarize(int nbonds,
472 const t_iatom forceatoms[],const t_iparams forceparams[],
473 const rvec x[],rvec f[],rvec fshift[],
474 const t_pbc *pbc,const t_graph *g,
475 real lambda,real *dvdlambda,
476 const t_mdatoms *md,t_fcdata *fcd,
477 int *global_atom_index)
478 {
479 int i,m,ki,ai,aj,type;
480 real dr,dr2,fbond,vbond,fij,vtot,ksh;
481 rvec dx;
482 ivec dt;
483
484 vtot = 0.0;
485 for(i=0; (i<nbonds); ) {
486 type = forceatoms[i++];
487 ai = forceatoms[i++];
488 aj = forceatoms[i++];
489 ksh = sqr(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;
490 if (debug)
491 fprintf(debug,"POL: local ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh);
492
493 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
494 dr2 = iprod(dx,dx); /* 5 */
495 dr = dr2*gmx_invsqrt(dr2); /* 10 */
496
497 *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond); /* 19 */
498
499 if (dr2 == 0.0)
500 continue;
501
502 vtot += vbond;/* 1*/
503 fbond *= gmx_invsqrt(dr2); /* 6 */
504
505 if (g) {
506 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
507 ki=IVEC2IS(dt);
508 }
509 for (m=0; (m<DIM); m++) { /* 15 */
510 fij=fbond*dx[m];
511 f[ai][m]+=fij;
512 f[aj][m]-=fij;
513 fshift[ki][m]+=fij;
514 fshift[CENTRAL][m]-=fij;
515 }
516 } /* 59 TOTAL */
517 return vtot;
518 }
519
520 real water_pol(int nbonds,
521 const t_iatom forceatoms[],const t_iparams forceparams[],
522 const rvec x[],rvec f[],rvec fshift[],
523 const t_pbc *pbc,const t_graph *g,
524 real lambda,real *dvdlambda,
525 const t_mdatoms *md,t_fcdata *fcd,
526 int *global_atom_index)
527 {
528 /* This routine implements anisotropic polarizibility for water, through
529 * a shell connected to a dummy with spring constant that differ in the
530 * three spatial dimensions in the molecular frame.
531 */
532 int i,m,aO,aH1,aH2,aD,aS,type,type0;
533 rvec dOH1,dOH2,dHH,dOD,dDS,nW,kk,dx,kdx,proj;
534 #ifdef DEBUG
535 rvec df;
536 #endif
537 real vtot,fij,r_HH,r_OD,r_nW,tx,ty,tz,qS;
538
539 vtot = 0.0;
540 if (nbonds > 0) {
541 type0 = forceatoms[0];
542 aS = forceatoms[5];
543 qS = md->chargeA[aS];
544 kk[XX] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_x;
545 kk[YY] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y;
546 kk[ZZ] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z;
547 r_HH = 1.0/forceparams[type0].wpol.rHH;
548 r_OD = 1.0/forceparams[type0].wpol.rOD;
549 if (debug) {
550 fprintf(debug,"WPOL: qS = %10.5f aS = %5d\n",qS,aS);
551 fprintf(debug,"WPOL: kk = %10.3f %10.3f %10.3f\n",
552 kk[XX],kk[YY],kk[ZZ]);
553 fprintf(debug,"WPOL: rOH = %10.3f rHH = %10.3f rOD = %10.3f\n",
554 forceparams[type0].wpol.rOH,
555 forceparams[type0].wpol.rHH,
556 forceparams[type0].wpol.rOD);
557 }
558 for(i=0; (i<nbonds); i+=6) {
559 type = forceatoms[i];
560 if (type != type0)
561 gmx_fatal(FARGS,"Sorry, type = %d, type0 = %d, file = %s, line = %d",
562 type,type0,__FILE__,__LINE__);
563 aO = forceatoms[i+1];
564 aH1 = forceatoms[i+2];
565 aH2 = forceatoms[i+3];
566 aD = forceatoms[i+4];
567 aS = forceatoms[i+5];
568
569 /* Compute vectors describing the water frame */
570 rvec_sub(x[aH1],x[aO], dOH1);
571 rvec_sub(x[aH2],x[aO], dOH2);
572 rvec_sub(x[aH2],x[aH1],dHH);
573 rvec_sub(x[aD], x[aO], dOD);
574 rvec_sub(x[aS], x[aD], dDS);
575 cprod(dOH1,dOH2,nW);
576
577 /* Compute inverse length of normal vector
578 * (this one could be precomputed, but I'm too lazy now)
579 */
580 r_nW = gmx_invsqrt(iprod(nW,nW));
581 /* This is for precision, but does not make a big difference,
582 * it can go later.
583 */
584 r_OD = gmx_invsqrt(iprod(dOD,dOD));
585
586 /* Normalize the vectors in the water frame */
587 svmul(r_nW,nW,nW);
588 svmul(r_HH,dHH,dHH);
589 svmul(r_OD,dOD,dOD);
590
591 /* Compute displacement of shell along components of the vector */
592 dx[ZZ] = iprod(dDS,dOD);
593 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
594 for(m=0; (m<DIM); m++)
595 proj[m] = dDS[m]-dx[ZZ]*dOD[m];
596
597 /*dx[XX] = iprod(dDS,nW);
598 dx[YY] = iprod(dDS,dHH);*/
599 dx[XX] = iprod(proj,nW);
600 for(m=0; (m<DIM); m++)
601 proj[m] -= dx[XX]*nW[m];
602 dx[YY] = iprod(proj,dHH);
603 /*#define DEBUG*/
604 #ifdef DEBUG
605 if (debug) {
606 fprintf(debug,"WPOL: dx2=%10g dy2=%10g dz2=%10g sum=%10g dDS^2=%10g\n",
607 sqr(dx[XX]),sqr(dx[YY]),sqr(dx[ZZ]),iprod(dx,dx),iprod(dDS,dDS));
608 fprintf(debug,"WPOL: dHH=(%10g,%10g,%10g)\n",dHH[XX],dHH[YY],dHH[ZZ]);
609 fprintf(debug,"WPOL: dOD=(%10g,%10g,%10g), 1/r_OD = %10g\n",
610 dOD[XX],dOD[YY],dOD[ZZ],1/r_OD);
611 fprintf(debug,"WPOL: nW =(%10g,%10g,%10g), 1/r_nW = %10g\n",
612 nW[XX],nW[YY],nW[ZZ],1/r_nW);
613 fprintf(debug,"WPOL: dx =%10g, dy =%10g, dz =%10g\n",
614 dx[XX],dx[YY],dx[ZZ]);
615 fprintf(debug,"WPOL: dDSx=%10g, dDSy=%10g, dDSz=%10g\n",
616 dDS[XX],dDS[YY],dDS[ZZ]);
617 }
618 #endif
619 /* Now compute the forces and energy */
620 kdx[XX] = kk[XX]*dx[XX];
621 kdx[YY] = kk[YY]*dx[YY];
622 kdx[ZZ] = kk[ZZ]*dx[ZZ];
623 vtot += iprod(dx,kdx);
624 for(m=0; (m<DIM); m++) {
625 /* This is a tensor operation but written out for speed */
626 tx = nW[m]*kdx[XX];
627 ty = dHH[m]*kdx[YY];
628 tz = dOD[m]*kdx[ZZ];
629 fij = -tx-ty-tz;
630 #ifdef DEBUG
631 df[m] = fij;
632 #endif
633 f[aS][m] += fij;
634 f[aD][m] -= fij;
635 }
636 #ifdef DEBUG
637 if (debug) {
638 fprintf(debug,"WPOL: vwpol=%g\n",0.5*iprod(dx,kdx));
639 fprintf(debug,"WPOL: df = (%10g, %10g, %10g)\n",df[XX],df[YY],df[ZZ]);
640 }
641 #endif
642 }
643 }
644 return 0.5*vtot;
645 }
646
647 static real do_1_thole(const rvec xi,const rvec xj,rvec fi,rvec fj,
648 const t_pbc *pbc,real qq,
649 rvec fshift[],real afac)
650 {
651 rvec r12;
652 real r12sq,r12_1,r12n,r12bar,v0,v1,fscal,ebar,fff;
653 int m,t;
654
655 t = pbc_rvec_sub(pbc,xi,xj,r12); /* 3 */
656
657 r12sq = iprod(r12,r12); /* 5 */
658 r12_1 = gmx_invsqrt(r12sq); /* 5 */
659 r12bar = afac/r12_1; /* 5 */
660 v0 = qq*ONE_4PI_EPS0*r12_1; /* 2 */
661 ebar = exp(-r12bar); /* 5 */
662 v1 = (1-(1+0.5*r12bar)*ebar); /* 4 */
663 fscal = ((v0*r12_1)*v1 - v0*0.5*afac*ebar*(r12bar+1))*r12_1; /* 9 */
664 if (debug)
665 fprintf(debug,"THOLE: v0 = %.3f v1 = %.3f r12= % .3f r12bar = %.3f fscal = %.3f ebar = %.3f\n",v0,v1,1/r12_1,r12bar,fscal,ebar);
666
667 for(m=0; (m<DIM); m++) {
668 fff = fscal*r12[m];
669 fi[m] += fff;
670 fj[m] -= fff;
671 fshift[t][m] += fff;
672 fshift[CENTRAL][m] -= fff;
673 } /* 15 */
674
675 return v0*v1; /* 1 */
676 /* 54 */
677 }
678
679 real thole_pol(int nbonds,
680 const t_iatom forceatoms[],const t_iparams forceparams[],
681 const rvec x[],rvec f[],rvec fshift[],
682 const t_pbc *pbc,const t_graph *g,
683 real lambda,real *dvdlambda,
684 const t_mdatoms *md,t_fcdata *fcd,
685 int *global_atom_index)
686 {
687 /* Interaction between two pairs of particles with opposite charge */
688 int i,type,a1,da1,a2,da2;
689 real q1,q2,qq,a,al1,al2,afac;
690 real V=0;
691
692 for(i=0; (i<nbonds); ) {
693 type = forceatoms[i++];
694 a1 = forceatoms[i++];
695 da1 = forceatoms[i++];
696 a2 = forceatoms[i++];
697 da2 = forceatoms[i++];
698 q1 = md->chargeA[da1];
699 q2 = md->chargeA[da2];
700 a = forceparams[type].thole.a;
701 al1 = forceparams[type].thole.alpha1;
702 al2 = forceparams[type].thole.alpha2;
703 qq = q1*q2;
704 afac = a*pow(al1*al2,-1.0/6.0);
705 V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq,fshift,afac);
706 V += do_1_thole(x[da1],x[a2], f[da1],f[a2], pbc,-qq,fshift,afac);
707 V += do_1_thole(x[a1], x[da2],f[a1], f[da2],pbc,-qq,fshift,afac);
708 V += do_1_thole(x[da1],x[da2],f[da1],f[da2],pbc, qq,fshift,afac);
709 }
710 /* 290 flops */
711 return V;
712 }
713
714 real bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
715 rvec r_ij,rvec r_kj,real *costh,
716 int *t1,int *t2)
717 /* Return value is the angle between the bonds i-j and j-k */
718 {
719 /* 41 FLOPS */
720 real th;
721
722 *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij); /* 3 */
723 *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj); /* 3 */
724
725 *costh=cos_angle(r_ij,r_kj); /* 25 */
726 th=acos(*costh); /* 10 */
727 /* 41 TOTAL */
728 return th;
729 }
730
731 real angles(int nbonds,
732 const t_iatom forceatoms[],const t_iparams forceparams[],
733 const rvec x[],rvec f[],rvec fshift[],
734 const t_pbc *pbc,const t_graph *g,
735 real lambda,real *dvdlambda,
736 const t_mdatoms *md,t_fcdata *fcd,
737 int *global_atom_index)
738 {
739 int i,ai,aj,ak,t1,t2,type;
740 rvec r_ij,r_kj;
741 real cos_theta,cos_theta2,theta,dVdt,va,vtot;
742 ivec jt,dt_ij,dt_kj;
743
744 vtot = 0.0;
745 for(i=0; (i<nbonds); ) {
746 type = forceatoms[i++];
747 ai = forceatoms[i++];
748 aj = forceatoms[i++];
749 ak = forceatoms[i++];
750
751 theta = bond_angle(x[ai],x[aj],x[ak],pbc,
752 r_ij,r_kj,&cos_theta,&t1,&t2); /* 41 */
753
754 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
755 forceparams[type].harmonic.krB,
756 forceparams[type].harmonic.rA*DEG2RAD,
757 forceparams[type].harmonic.rB*DEG2RAD,
758 theta,lambda,&va,&dVdt); /* 21 */
759 vtot += va;
760
761 cos_theta2 = sqr(cos_theta);
762 if (cos_theta2 < 1) {
763 int m;
764 real st,sth;
765 real cik,cii,ckk;
766 real nrkj2,nrij2;
767 rvec f_i,f_j,f_k;
768
769 st = dVdt*gmx_invsqrt(1 - cos_theta2); /* 12 */
770 sth = st*cos_theta; /* 1 */
771 #ifdef DEBUG
772 if (debug)
773 fprintf(debug,"ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n",
774 theta*RAD2DEG,va,dVdt);
775 #endif
776 nrkj2=iprod(r_kj,r_kj); /* 5 */
777 nrij2=iprod(r_ij,r_ij);
778
779 cik=st*gmx_invsqrt(nrkj2*nrij2); /* 12 */
780 cii=sth/nrij2; /* 10 */
781 ckk=sth/nrkj2; /* 10 */
782
783 for (m=0; (m<DIM); m++) { /* 39 */
784 f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
785 f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
786 f_j[m]=-f_i[m]-f_k[m];
787 f[ai][m]+=f_i[m];
788 f[aj][m]+=f_j[m];
789 f[ak][m]+=f_k[m];
790 }
791 if (g) {
792 copy_ivec(SHIFT_IVEC(g,aj),jt);
793
794 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
795 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
796 t1=IVEC2IS(dt_ij);
797 t2=IVEC2IS(dt_kj);
798 }
799 rvec_inc(fshift[t1],f_i);
800 rvec_inc(fshift[CENTRAL],f_j);
801 rvec_inc(fshift[t2],f_k);
802 } /* 161 TOTAL */
803 }
804 return vtot;
805 }
806
807 real urey_bradley(int nbonds,
808 const t_iatom forceatoms[],const t_iparams forceparams[],
809 const rvec x[],rvec f[],rvec fshift[],
810 const t_pbc *pbc,const t_graph *g,
811 real lambda,real *dvdlambda,
812 const t_mdatoms *md,t_fcdata *fcd,
813 int *global_atom_index)
814 {
815 int i,m,ai,aj,ak,t1,t2,type,ki;
816 rvec r_ij,r_kj,r_ik;
817 real cos_theta,cos_theta2,theta;
818 real dVdt,va,vtot,kth,th0,kUB,r13,dr,dr2,vbond,fbond,fik;
819 ivec jt,dt_ij,dt_kj,dt_ik;
820
821 vtot = 0.0;
822 for(i=0; (i<nbonds); ) {
823 type = forceatoms[i++];
824 ai = forceatoms[i++];
825 aj = forceatoms[i++];
826 ak = forceatoms[i++];
827 th0 = forceparams[type].u_b.theta*DEG2RAD;
828 kth = forceparams[type].u_b.ktheta;
829 r13 = forceparams[type].u_b.r13;
830 kUB = forceparams[type].u_b.kUB;
831
832 theta = bond_angle(x[ai],x[aj],x[ak],pbc,
833 r_ij,r_kj,&cos_theta,&t1,&t2); /* 41 */
834
835 *dvdlambda += harmonic(kth,kth,th0,th0,theta,lambda,&va,&dVdt); /* 21 */
836 vtot += va;
837
838 ki = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik); /* 3 */
839 dr2 = iprod(r_ik,r_ik); /* 5 */
840 dr = dr2*gmx_invsqrt(dr2); /* 10 */
841
842 *dvdlambda += harmonic(kUB,kUB,r13,r13,dr,lambda,&vbond,&fbond); /* 19 */
843
844 cos_theta2 = sqr(cos_theta); /* 1 */
845 if (cos_theta2 < 1) {
846 real st,sth;
847 real cik,cii,ckk;
848 real nrkj2,nrij2;
849 rvec f_i,f_j,f_k;
850
851 st = dVdt*gmx_invsqrt(1 - cos_theta2); /* 12 */
852 sth = st*cos_theta; /* 1 */
853 #ifdef DEBUG
854 if (debug)
855 fprintf(debug,"ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n",
856 theta*RAD2DEG,va,dVdt);
857 #endif
858 nrkj2=iprod(r_kj,r_kj); /* 5 */
859 nrij2=iprod(r_ij,r_ij);
860
861 cik=st*gmx_invsqrt(nrkj2*nrij2); /* 12 */
862 cii=sth/nrij2; /* 10 */
863 ckk=sth/nrkj2; /* 10 */
864
865 for (m=0; (m<DIM); m++) { /* 39 */
866 f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
867 f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
868 f_j[m]=-f_i[m]-f_k[m];
869 f[ai][m]+=f_i[m];
870 f[aj][m]+=f_j[m];
871 f[ak][m]+=f_k[m];
872 }
873 if (g) {
874 copy_ivec(SHIFT_IVEC(g,aj),jt);
875
876 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
877 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
878 t1=IVEC2IS(dt_ij);
879 t2=IVEC2IS(dt_kj);
880 }
881 rvec_inc(fshift[t1],f_i);
882 rvec_inc(fshift[CENTRAL],f_j);
883 rvec_inc(fshift[t2],f_k);
884 } /* 161 TOTAL */
885 /* Time for the bond calculations */
886 if (dr2 == 0.0)
887 continue;
888
889 vtot += vbond; /* 1*/
890 fbond *= gmx_invsqrt(dr2); /* 6 */
891
892 if (g) {
893 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),dt_ik);
894 ki=IVEC2IS(dt_ik);
895 }
896 for (m=0; (m<DIM); m++) { /* 15 */
897 fik=fbond*r_ik[m];
898 f[ai][m]+=fik;
899 f[ak][m]-=fik;
900 fshift[ki][m]+=fik;
901 fshift[CENTRAL][m]-=fik;
902 }
903 }
904 return vtot;
905 }
906
907 real quartic_angles(int nbonds,
908 const t_iatom forceatoms[],const t_iparams forceparams[],
909 const rvec x[],rvec f[],rvec fshift[],
910 const t_pbc *pbc,const t_graph *g,
911 real lambda,real *dvdlambda,
912 const t_mdatoms *md,t_fcdata *fcd,
913 int *global_atom_index)
914 {
915 int i,j,ai,aj,ak,t1,t2,type;
916 rvec r_ij,r_kj;
917 real cos_theta,cos_theta2,theta,dt,dVdt,va,dtp,c,vtot;
918 ivec jt,dt_ij,dt_kj;
919
920 vtot = 0.0;
921 for(i=0; (i<nbonds); ) {
922 type = forceatoms[i++];
923 ai = forceatoms[i++];
924 aj = forceatoms[i++];
925 ak = forceatoms[i++];
926
927 theta = bond_angle(x[ai],x[aj],x[ak],pbc,
928 r_ij,r_kj,&cos_theta,&t1,&t2); /* 41 */
929
930 dt = theta - forceparams[type].qangle.theta*DEG2RAD; /* 2 */
931
932 dVdt = 0;
933 va = forceparams[type].qangle.c[0];
934 dtp = 1.0;
935 for(j=1; j<=4; j++) {
936 c = forceparams[type].qangle.c[j];
937 dVdt -= j*c*dtp;
938 dtp *= dt;
939 va += c*dtp;
940 }
941 /* 20 */
942
943 vtot += va;
944
945 cos_theta2 = sqr(cos_theta); /* 1 */
946 if (cos_theta2 < 1) {
947 int m;
948 real st,sth;
949 real cik,cii,ckk;
950 real nrkj2,nrij2;
951 rvec f_i,f_j,f_k;
952
953 st = dVdt*gmx_invsqrt(1 - cos_theta2); /* 12 */
954 sth = st*cos_theta; /* 1 */
955 #ifdef DEBUG
956 if (debug)
957 fprintf(debug,"ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n",
958 theta*RAD2DEG,va,dVdt);
959 #endif
960 nrkj2=iprod(r_kj,r_kj); /* 5 */
961 nrij2=iprod(r_ij,r_ij);
962
963 cik=st*gmx_invsqrt(nrkj2*nrij2); /* 12 */
964 cii=sth/nrij2; /* 10 */
965 ckk=sth/nrkj2; /* 10 */
966
967 for (m=0; (m<DIM); m++) { /* 39 */
968 f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
969 f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
970 f_j[m]=-f_i[m]-f_k[m];
971 f[ai][m]+=f_i[m];
972 f[aj][m]+=f_j[m];
973 f[ak][m]+=f_k[m];
974 }
975 if (g) {
976 copy_ivec(SHIFT_IVEC(g,aj),jt);
977
978 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
979 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
980 t1=IVEC2IS(dt_ij);
981 t2=IVEC2IS(dt_kj);
982 }
983 rvec_inc(fshift[t1],f_i);
984 rvec_inc(fshift[CENTRAL],f_j);
985 rvec_inc(fshift[t2],f_k);
986 } /* 153 TOTAL */
987 }
988 return vtot;
989 }
990
991 real dih_angle(const rvec xi,const rvec xj,const rvec xk,const rvec xl,
992 const t_pbc *pbc,
993 rvec r_ij,rvec r_kj,rvec r_kl,rvec m,rvec n,
994 real *sign,int *t1,int *t2,int *t3)
995 {
996 real ipr,phi;
997
998 *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij); /* 3 */
999 *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj); /* 3 */
1000 *t3 = pbc_rvec_sub(pbc,xk,xl,r_kl); /* 3 */
1001
1002 cprod(r_ij,r_kj,m); /* 9 */
1003 cprod(r_kj,r_kl,n); /* 9 */
1004 phi=gmx_angle(m,n); /* 49 (assuming 25 for atan2) */
1005 ipr=iprod(r_ij,n); /* 5 */
1006 (*sign)=(ipr<0.0)?-1.0:1.0;
1007 phi=(*sign)*phi; /* 1 */
1008 /* 82 TOTAL */
1009 return phi;
1010 }
1011
1012
1013
1014 void do_dih_fup(int i,int j,int k,int l,real ddphi,
1015 rvec r_ij,rvec r_kj,rvec r_kl,
1016 rvec m,rvec n,rvec f[],rvec fshift[],
1017 const t_pbc *pbc,const t_graph *g,
1018 const rvec x[],int t1,int t2,int t3)
1019 {
1020 /* 143 FLOPS */
1021 rvec f_i,f_j,f_k,f_l;
1022 rvec uvec,vvec,svec,dx_jl;
1023 real iprm,iprn,nrkj,nrkj2;
1024 real a,p,q,toler;
1025 ivec jt,dt_ij,dt_kj,dt_lj;
1026
1027 iprm = iprod(m,m); /* 5 */
1028 iprn = iprod(n,n); /* 5 */
1029 nrkj2 = iprod(r_kj,r_kj); /* 5 */
1030 toler = nrkj2*GMX_REAL_EPS;
1031 if ((iprm > toler) && (iprn > toler)) {
1032 nrkj = nrkj2*gmx_invsqrt(nrkj2); /* 10 */
1033 a = -ddphi*nrkj/iprm; /* 11 */
1034 svmul(a,m,f_i); /* 3 */
1035 a = ddphi*nrkj/iprn; /* 11 */
1036 svmul(a,n,f_l); /* 3 */
1037 p = iprod(r_ij,r_kj); /* 5 */
1038 p /= nrkj2; /* 10 */
1039 q = iprod(r_kl,r_kj); /* 5 */
1040 q /= nrkj2; /* 10 */
1041 svmul(p,f_i,uvec); /* 3 */
1042 svmul(q,f_l,vvec); /* 3 */
1043 rvec_sub(uvec,vvec,svec); /* 3 */
1044 rvec_sub(f_i,svec,f_j); /* 3 */
1045 rvec_add(f_l,svec,f_k); /* 3 */
1046 rvec_inc(f[i],f_i); /* 3 */
1047 rvec_dec(f[j],f_j); /* 3 */
1048 rvec_dec(f[k],f_k); /* 3 */
1049 rvec_inc(f[l],f_l); /* 3 */
1050
1051 if (g) {
1052 copy_ivec(SHIFT_IVEC(g,j),jt);
1053 ivec_sub(SHIFT_IVEC(g,i),jt,dt_ij);
1054 ivec_sub(SHIFT_IVEC(g,k),jt,dt_kj);
1055 ivec_sub(SHIFT_IVEC(g,l),jt,dt_lj);
1056 t1=IVEC2IS(dt_ij);
1057 t2=IVEC2IS(dt_kj);
1058 t3=IVEC2IS(dt_lj);
1059 } else if (pbc) {
1060 t3 = pbc_rvec_sub(pbc,x[l],x[j],dx_jl);
1061 } else {
1062 t3 = CENTRAL;
1063 }
1064
1065 rvec_inc(fshift[t1],f_i);
1066 rvec_dec(fshift[CENTRAL],f_j);
1067 rvec_dec(fshift[t2],f_k);
1068 rvec_inc(fshift[t3],f_l);
1069 }
1070 /* 112 TOTAL */
1071 }
1072
1073
1074 real dopdihs(real cpA,real cpB,real phiA,real phiB,int mult,
1075 real phi,real lambda,real *V,real *F)
1076 {
1077 real v,dvdl,mdphi,v1,sdphi,ddphi;
1078 real L1 = 1.0 - lambda;
1079 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1080 real dph0 = (phiB - phiA)*DEG2RAD;
1081 real cp = L1*cpA + lambda*cpB;
1082
1083 mdphi = mult*phi - ph0;
1084 sdphi = sin(mdphi);
1085 ddphi = -cp*mult*sdphi;
1086 v1 = 1.0 + cos(mdphi);
1087 v = cp*v1;
1088
1089 dvdl = (cpB - cpA)*v1 + cp*dph0*sdphi;
1090
1091 *V = v;
1092 *F = ddphi;
1093
1094 return dvdl;
1095
1096 /* That was 40 flops */
1097 }
1098
1099 static real dopdihs_min(real cpA,real cpB,real phiA,real phiB,int mult,
1100 real phi,real lambda,real *V,real *F)
1101 /* similar to dopdihs, except for a minus sign *
1102 * and a different treatment of mult/phi0 */
1103 {
1104 real v,dvdl,mdphi,v1,sdphi,ddphi;
1105 real L1 = 1.0 - lambda;
1106 real ph0 = (L1*phiA + lambda*phiB)*DEG2RAD;
1107 real dph0 = (phiB - phiA)*DEG2RAD;
1108 real cp = L1*cpA + lambda*cpB;
1109
1110 mdphi = mult*(phi-ph0);
1111 sdphi = sin(mdphi);
1112 ddphi = cp*mult*sdphi;
1113 v1 = 1.0-cos(mdphi);
1114 v = cp*v1;
1115
1116 dvdl = (cpB-cpA)*v1 + cp*dph0*sdphi;
1117
1118 *V = v;
1119 *F = ddphi;
1120
1121 return dvdl;
1122
1123 /* That was 40 flops */
1124 }
1125
1126 real pdihs(int nbonds,
1127 const t_iatom forceatoms[],const t_iparams forceparams[],
1128 const rvec x[],rvec f[],rvec fshift[],
1129 const t_pbc *pbc,const t_graph *g,
1130 real lambda,real *dvdlambda,
1131 const t_mdatoms *md,t_fcdata *fcd,
1132 int *global_atom_index)
1133 {
1134 int i,type,ai,aj,ak,al;
1135 int t1,t2,t3;
1136 rvec r_ij,r_kj,r_kl,m,n;
1137 real phi,sign,ddphi,vpd,vtot;
1138
1139 vtot = 0.0;
1140
1141 for(i=0; (i<nbonds); ) {
1142 type = forceatoms[i++];
1143 ai = forceatoms[i++];
1144 aj = forceatoms[i++];
1145 ak = forceatoms[i++];
1146 al = forceatoms[i++];
1147
1148 phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1149 &sign,&t1,&t2,&t3); /* 84 */
1150
1151 *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
1152 forceparams[type].pdihs.cpB,
1153 forceparams[type].pdihs.phiA,
1154 forceparams[type].pdihs.phiB,
1155 forceparams[type].pdihs.mult,
1156 phi,lambda,&vpd,&ddphi);
1157
1158 vtot += vpd;
1159 do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1160 f,fshift,pbc,g,x,t1,t2,t3); /* 112 */
1161
1162 #ifdef DEBUG
1163 fprintf(debug,"pdih: (%d,%d,%d,%d) phi=%g\n",
1164 ai,aj,ak,al,phi);
1165 #endif
1166 } /* 223 TOTAL */
1167
1168 return vtot;
1169 }
1170
1171
1172
1173 real idihs(int nbonds,
1174 const t_iatom forceatoms[],const t_iparams forceparams[],
1175 const rvec x[],rvec f[],rvec fshift[],
1176 const t_pbc *pbc,const t_graph *g,
1177 real lambda,real *dvdlambda,
1178 const t_mdatoms *md,t_fcdata *fcd,
1179 int *global_atom_index)
1180 {
1181 int i,type,ai,aj,ak,al;
1182 int t1,t2,t3;
1183 real phi,phi0,dphi0,ddphi,sign,vtot;
1184 rvec r_ij,r_kj,r_kl,m,n;
1185 real L1,kk,dp,dp2,kA,kB,pA,pB,dvdl;
1186
1187 L1 = 1.0-lambda;
1188 dvdl = 0;
1189
1190 vtot = 0.0;
1191 for(i=0; (i<nbonds); ) {
1192 type = forceatoms[i++];
1193 ai = forceatoms[i++];
1194 aj = forceatoms[i++];
1195 ak = forceatoms[i++];
1196 al = forceatoms[i++];
1197
1198 phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1199 &sign,&t1,&t2,&t3); /* 84 */
1200
1201 /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
1202 * force changes if we just apply a normal harmonic.
1203 * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
1204 * This means we will never have the periodicity problem, unless
1205 * the dihedral is Pi away from phiO, which is very unlikely due to
1206 * the potential.
1207 */
1208 kA = forceparams[type].harmonic.krA;
1209 kB = forceparams[type].harmonic.krB;
1210 pA = forceparams[type].harmonic.rA;
1211 pB = forceparams[type].harmonic.rB;
1212
1213 kk = L1*kA + lambda*kB;
1214 phi0 = (L1*pA + lambda*pB)*DEG2RAD;
1215 dphi0 = (pB - pA)*DEG2RAD;
1216
1217 /* dp = (phi-phi0), modulo (-pi,pi) */
1218 dp = phi-phi0;
1219 /* dp cannot be outside (-2*pi,2*pi) */
1220 if (dp >= M_PI)
1221 dp -= 2*M_PI;
1222 else if(dp < -M_PI)
1223 dp += 2*M_PI;
1224
1225 dp2 = dp*dp;
1226
1227 vtot += 0.5*kk*dp2;
1228 ddphi = -kk*dp;
1229
1230 dvdl += 0.5*(kB - kA)*dp2 - kk*dphi0*dp;
1231
1232 do_dih_fup(ai,aj,ak,al,(real)(-ddphi),r_ij,r_kj,r_kl,m,n,
1233 f,fshift,pbc,g,x,t1,t2,t3); /* 112 */
1234 /* 217 TOTAL */
1235 #ifdef DEBUG
1236 if (debug)
1237 fprintf(debug,"idih: (%d,%d,%d,%d) phi=%g\n",
1238 ai,aj,ak,al,phi);
1239 #endif
1240 }
1241
1242 *dvdlambda += dvdl;
1243 return vtot;
1244 }
1245
1246
1247 real posres(int nbonds,
1248 const t_iatom forceatoms[],const t_iparams forceparams[],
1249 const rvec x[],rvec f[],rvec vir_diag,
1250 t_pbc *pbc,
1251 real lambda,real *dvdlambda,
1252 int refcoord_scaling,int ePBC,rvec comA,rvec comB)
1253 {
1254 int i,ai,m,d,type,ki,npbcdim=0;
1255 const t_iparams *pr;
1256 real L1;
1257 real vtot,kk,fm;
1258 real posA,posB,ref=0;
1259 rvec comA_sc,comB_sc,rdist,dpdl,pos,dx;
1260
1261 npbcdim = ePBC2npbcdim(ePBC);
1262
1263 if (refcoord_scaling == erscCOM)
1264 {
1265 clear_rvec(comA_sc);
1266 clear_rvec(comB_sc);
1267 for(m=0; m<npbcdim; m++)
1268 {
1269 for(d=m; d<npbcdim; d++)
1270 {
1271 comA_sc[m] += comA[d]*pbc->box[d][m];
1272 comB_sc[m] += comB[d]*pbc->box[d][m];
1273 }
1274 }
1275 }
1276
1277 L1 = 1.0 - lambda;
1278
1279 vtot = 0.0;
1280 for(i=0; (i<nbonds); )
1281 {
1282 type = forceatoms[i++];
1283 ai = forceatoms[i++];
1284 pr = &forceparams[type];
1285
1286 for(m=0; m<DIM; m++)
1287 {
1288 posA = forceparams[type].posres.pos0A[m];
1289 posB = forceparams[type].posres.pos0B[m];
1290 if (m < npbcdim)
1291 {
1292 switch (refcoord_scaling)
1293 {
1294 case erscNO:
1295 ref = 0;
1296 rdist[m] = L1*posA + lambda*posB;
1297 dpdl[m] = posB - posA;
1298 break;
1299 case erscALL:
1300 /* Box relative coordinates are stored for dimensions with pbc */
1301 posA *= pbc->box[m][m];
1302 posB *= pbc->box[m][m];
1303 for(d=m+1; d<npbcdim; d++)
1304 {
1305 posA += forceparams[type].posres.pos0A[d]*pbc->box[d][m];
1306 posB += forceparams[type].posres.pos0B[d]*pbc->box[d][m];
1307 }
1308 ref = L1*posA + lambda*posB;
1309 rdist[m] = 0;
1310 dpdl[m] = posB - posA;
1311 break;
1312 case erscCOM:
1313 ref = L1*comA_sc[m] + lambda*comB_sc[m];
1314 rdist[m] = L1*posA + lambda*posB;
1315 dpdl[m] = comB_sc[m] - comA_sc[m] + posB - posA;
1316 break;
1317 }
1318 }
1319 else
1320 {
1321 ref = L1*posA + lambda*posB;
1322 rdist[m] = 0;
1323 dpdl[m] = posB - posA;
1324 }
1325
1326 /* We do pbc_dx with ref+rdist,
1327 * since with only ref we can be up to half a box vector wrong.
1328 */
1329 pos[m] = ref + rdist[m];
1330 }
1331
1332 if (pbc)
1333 {
1334 pbc_dx(pbc,x[ai],pos,dx);
1335 }
1336 else
1337 {
1338 rvec_sub(x[ai],pos,dx);
1339 }
1340
1341 for (m=0; (m<DIM); m++)
1342 {
1343 kk = L1*pr->posres.fcA[m] + lambda*pr->posres.fcB[m];
1344 fm = -kk*dx[m];
1345 f[ai][m] += fm;
1346 vtot += 0.5*kk*dx[m]*dx[m];
1347 *dvdlambda +=
1348 0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
1349 -fm*dpdl[m];
1350
1351 /* Here we correct for the pbc_dx which included rdist */
1352 vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm;
1353 }
1354 }
1355
1356 return vtot;
1357 }
1358
1359 static real low_angres(int nbonds,
1360 const t_iatom forceatoms[],const t_iparams forceparams[],
1361 const rvec x[],rvec f[],rvec fshift[],
1362 const t_pbc *pbc,const t_graph *g,
1363 real lambda,real *dvdlambda,
1364 gmx_bool bZAxis)
1365 {
1366 int i,m,type,ai,aj,ak,al;
1367 int t1,t2;
1368 real phi,cos_phi,cos_phi2,vid,vtot,dVdphi;
1369 rvec r_ij,r_kl,f_i,f_k={0,0,0};
1370 real st,sth,nrij2,nrkl2,c,cij,ckl;
1371
1372 ivec dt;
1373 t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
1374
1375 vtot = 0.0;
1376 ak=al=0; /* to avoid warnings */
1377 for(i=0; i<nbonds; ) {
1378 type = forceatoms[i++];
1379 ai = forceatoms[i++];
1380 aj = forceatoms[i++];
1381 t1 = pbc_rvec_sub(pbc,x[aj],x[ai],r_ij); /* 3 */
1382 if (!bZAxis) {
1383 ak = forceatoms[i++];
1384 al = forceatoms[i++];
1385 t2 = pbc_rvec_sub(pbc,x[al],x[ak],r_kl); /* 3 */
1386 } else {
1387 r_kl[XX] = 0;
1388 r_kl[YY] = 0;
1389 r_kl[ZZ] = 1;
1390 }
1391
1392 cos_phi = cos_angle(r_ij,r_kl); /* 25 */
1393 phi = acos(cos_phi); /* 10 */
1394
1395 *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
1396 forceparams[type].pdihs.cpB,
1397 forceparams[type].pdihs.phiA,
1398 forceparams[type].pdihs.phiB,
1399 forceparams[type].pdihs.mult,
1400 phi,lambda,&vid,&dVdphi); /* 40 */
1401
1402 vtot += vid;
1403
1404 cos_phi2 = sqr(cos_phi); /* 1 */
1405 if (cos_phi2 < 1) {
1406 st = -dVdphi*gmx_invsqrt(1 - cos_phi2); /* 12 */
1407 sth = st*cos_phi; /* 1 */
1408 nrij2 = iprod(r_ij,r_ij); /* 5 */
1409 nrkl2 = iprod(r_kl,r_kl); /* 5 */
1410
1411 c = st*gmx_invsqrt(nrij2*nrkl2); /* 11 */
1412 cij = sth/nrij2; /* 10 */
1413 ckl = sth/nrkl2; /* 10 */
1414
1415 for (m=0; m<DIM; m++) { /* 18+18 */
1416 f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
1417 f[ai][m] += f_i[m];
1418 f[aj][m] -= f_i[m];
1419 if (!bZAxis) {
1420 f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
1421 f[ak][m] += f_k[m];
1422 f[al][m] -= f_k[m];
1423 }
1424 }
1425
1426 if (g) {
1427 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1428 t1=IVEC2IS(dt);
1429 }
1430 rvec_inc(fshift[t1],f_i);
1431 rvec_dec(fshift[CENTRAL],f_i);
1432 if (!bZAxis) {
1433 if (g) {
1434 ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,al),dt);
1435 t2=IVEC2IS(dt);
1436 }
1437 rvec_inc(fshift[t2],f_k);
1438 rvec_dec(fshift[CENTRAL],f_k);
1439 }
1440 }
1441 }
1442
1443 return vtot; /* 184 / 157 (bZAxis) total */
1444 }
1445
1446 real angres(int nbonds,
1447 const t_iatom forceatoms[],const t_iparams forceparams[],
1448 const rvec x[],rvec f[],rvec fshift[],
1449 const t_pbc *pbc,const t_graph *g,
1450 real lambda,real *dvdlambda,
1451 const t_mdatoms *md,t_fcdata *fcd,
1452 int *global_atom_index)
1453 {
1454 return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1455 lambda,dvdlambda,FALSE);
1456 }
1457
1458 real angresz(int nbonds,
1459 const t_iatom forceatoms[],const t_iparams forceparams[],
1460 const rvec x[],rvec f[],rvec fshift[],
1461 const t_pbc *pbc,const t_graph *g,
1462 real lambda,real *dvdlambda,
1463 const t_mdatoms *md,t_fcdata *fcd,
1464 int *global_atom_index)
1465 {
1466 return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1467 lambda,dvdlambda,TRUE);
1468 }
1469
1470
1471 real unimplemented(int nbonds,
1472 const t_iatom forceatoms[],const t_iparams forceparams[],
1473 const rvec x[],rvec f[],rvec fshift[],
1474 const t_pbc *pbc,const t_graph *g,
1475 real lambda,real *dvdlambda,
1476 const t_mdatoms *md,t_fcdata *fcd,
1477 int *global_atom_index)
1478 {
1479 gmx_impl("*** you are using a not implemented function");
1480
1481 return 0.0; /* To make the compiler happy */
1482 }
1483
1484 real rbdihs(int nbonds,
1485 const t_iatom forceatoms[],const t_iparams forceparams[],
1486 const rvec x[],rvec f[],rvec fshift[],
1487 const t_pbc *pbc,const t_graph *g,
1488 real lambda,real *dvdlambda,
1489 const t_mdatoms *md,t_fcdata *fcd,
1490 int *global_atom_index)
1491 {
1492 const real c0=0.0,c1=1.0,c2=2.0,c3=3.0,c4=4.0,c5=5.0;
1493 int type,ai,aj,ak,al,i,j;
1494 int t1,t2,t3;
1495 rvec r_ij,r_kj,r_kl,m,n;
1496 real parmA[NR_RBDIHS];
1497 real parmB[NR_RBDIHS];
1498 real parm[NR_RBDIHS];
1499 real cos_phi,phi,rbp,rbpBA;
1500 real v,sign,ddphi,sin_phi;
1501 real cosfac,vtot;
1502 real L1 = 1.0-lambda;
1503 real dvdl=0;
1504
1505 vtot = 0.0;
1506 for(i=0; (i<nbonds); ) {
1507 type = forceatoms[i++];
1508 ai = forceatoms[i++];
1509 aj = forceatoms[i++];
1510 ak = forceatoms[i++];
1511 al = forceatoms[i++];
1512
1513 phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1514 &sign,&t1,&t2,&t3); /* 84 */
1515
1516 /* Change to polymer convention */
1517 if (phi < c0)
1518 phi += M_PI;
1519 else
1520 phi -= M_PI; /* 1 */
1521
1522 cos_phi = cos(phi);
1523 /* Beware of accuracy loss, cannot use 1-sqrt(cos^2) ! */
1524 sin_phi = sin(phi);
1525
1526 for(j=0; (j<NR_RBDIHS); j++) {
1527 parmA[j] = forceparams[type].rbdihs.rbcA[j];
1528 parmB[j] = forceparams[type].rbdihs.rbcB[j];
1529 parm[j] = L1*parmA[j]+lambda*parmB[j];
1530 }
1531 /* Calculate cosine powers */
1532 /* Calculate the energy */
1533 /* Calculate the derivative */
1534
1535 v = parm[0];
1536 dvdl += (parmB[0]-parmA[0]);
1537 ddphi = c0;
1538 cosfac = c1;
1539
1540 rbp = parm[1];
1541 rbpBA = parmB[1]-parmA[1];
1542 ddphi += rbp*cosfac;
1543 cosfac *= cos_phi;
1544 v += cosfac*rbp;
1545 dvdl += cosfac*rbpBA;
1546 rbp = parm[2];
1547 rbpBA = parmB[2]-parmA[2];
1548 ddphi += c2*rbp*cosfac;
1549 cosfac *= cos_phi;
1550 v += cosfac*rbp;
1551 dvdl += cosfac*rbpBA;
1552 rbp = parm[3];
1553 rbpBA = parmB[3]-parmA[3];
1554 ddphi += c3*rbp*cosfac;
1555 cosfac *= cos_phi;
1556 v += cosfac*rbp;
1557 dvdl += cosfac*rbpBA;
1558 rbp = parm[4];
1559 rbpBA = parmB[4]-parmA[4];
1560 ddphi += c4*rbp*cosfac;
1561 cosfac *= cos_phi;
1562 v += cosfac*rbp;
1563 dvdl += cosfac*rbpBA;
1564 rbp = parm[5];
1565 rbpBA = parmB[5]-parmA[5];
1566 ddphi += c5*rbp*cosfac;
1567 cosfac *= cos_phi;
1568 v += cosfac*rbp;
1569 dvdl += cosfac*rbpBA;
1570
1571 ddphi = -ddphi*sin_phi; /* 11 */
1572
1573 do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1574 f,fshift,pbc,g,x,t1,t2,t3); /* 112 */
1575 vtot += v;
1576 }
1577 *dvdlambda += dvdl;
1578
1579 return vtot;
1580 }
1581
1582 int cmap_setup_grid_index(int ip, int grid_spacing, int *ipm1, int *ipp1, int *ipp2)
1583 {
1584 int im1, ip1, ip2;
1585
1586 if(ip<0)
1587 {
1588 ip = ip + grid_spacing - 1;
1589 }
1590 else if(ip > grid_spacing)
1591 {
1592 ip = ip - grid_spacing - 1;
1593 }
1594
1595 im1 = ip - 1;
1596 ip1 = ip + 1;
1597 ip2 = ip + 2;
1598
1599 if(ip == 0)
1600 {
1601 im1 = grid_spacing - 1;
1602 }
1603 else if(ip == grid_spacing-2)
1604 {
1605 ip2 = 0;
1606 }
1607 else if(ip == grid_spacing-1)
1608 {
1609 ip1 = 0;
1610 ip2 = 1;
1611 }
1612
1613 *ipm1 = im1;
1614 *ipp1 = ip1;
1615 *ipp2 = ip2;
1616
1617 return ip;
1618
1619 }
1620
1621 real cmap_dihs(int nbonds,
1622 const t_iatom forceatoms[],const t_iparams forceparams[],
1623 const gmx_cmap_t *cmap_grid,
1624 const rvec x[],rvec f[],rvec fshift[],
1625 const t_pbc *pbc,const t_graph *g,
1626 real lambda,real *dvdlambda,
1627 const t_mdatoms *md,t_fcdata *fcd,
1628 int *global_atom_index)
1629 {
1630 int i,j,k,n,idx;
1631 int ai,aj,ak,al,am;
1632 int a1i,a1j,a1k,a1l,a2i,a2j,a2k,a2l;
1633 int type,cmapA;
1634 int t11,t21,t31,t12,t22,t32;
1635 int iphi1,ip1m1,ip1p1,ip1p2;
1636 int iphi2,ip2m1,ip2p1,ip2p2;
1637 int l1,l2,l3,l4;
1638 int pos1,pos2,pos3,pos4,tmp;
1639
1640 real ty[4],ty1[4],ty2[4],ty12[4],tc[16],tx[16];
1641 real phi1,psi1,cos_phi1,sin_phi1,sign1,xphi1;
1642 real phi2,psi2,cos_phi2,sin_phi2,sign2,xphi2;
1643 real dx,xx,tt,tu,e,df1,df2,ddf1,ddf2,ddf12,vtot;
1644 real ra21,rb21,rg21,rg1,rgr1,ra2r1,rb2r1,rabr1;
1645 real ra22,rb22,rg22,rg2,rgr2,ra2r2,rb2r2,rabr2;
1646 real fg1,hg1,fga1,hgb1,gaa1,gbb1;
1647 real fg2,hg2,fga2,hgb2,gaa2,gbb2;
1648 real fac;
1649
1650 rvec r1_ij, r1_kj, r1_kl,m1,n1;
1651 rvec r2_ij, r2_kj, r2_kl,m2,n2;
1652 rvec f1_i,f1_j,f1_k,f1_l;
1653 rvec f2_i,f2_j,f2_k,f2_l;
1654 rvec a1,b1,a2,b2;
1655 rvec f1,g1,h1,f2,g2,h2;
1656 rvec dtf1,dtg1,dth1,dtf2,dtg2,dth2;
1657 ivec jt1,dt1_ij,dt1_kj,dt1_lj;
1658 ivec jt2,dt2_ij,dt2_kj,dt2_lj;
1659
1660 const real *cmapd;
1661
1662 int loop_index[4][4] = {
1663 {0,4,8,12},
1664 {1,5,9,13},
1665 {2,6,10,14},
1666 {3,7,11,15}
1667 };
1668
1669 /* Total CMAP energy */
1670 vtot = 0;
1671
1672 for(n=0;n<nbonds; )
1673 {
1674 /* Five atoms are involved in the two torsions */
1675 type = forceatoms[n++];
1676 ai = forceatoms[n++];
1677 aj = forceatoms[n++];
1678 ak = forceatoms[n++];
1679 al = forceatoms[n++];
1680 am = forceatoms[n++];
1681
1682 /* Which CMAP type is this */
1683 cmapA = forceparams[type].cmap.cmapA;
1684 cmapd = cmap_grid->cmapdata[cmapA].cmap;
1685
1686 /* First torsion */
1687 a1i = ai;
1688 a1j = aj;
1689 a1k = ak;
1690 a1l = al;
1691
1692 phi1 = dih_angle(x[a1i], x[a1j], x[a1k], x[a1l], pbc, r1_ij, r1_kj, r1_kl, m1, n1,
1693 &sign1, &t11, &t21, &t31); /* 84 */
1694
1695 cos_phi1 = cos(phi1);
1696
1697 a1[0] = r1_ij[1]*r1_kj[2]-r1_ij[2]*r1_kj[1];
1698 a1[1] = r1_ij[2]*r1_kj[0]-r1_ij[0]*r1_kj[2];
1699 a1[2] = r1_ij[0]*r1_kj[1]-r1_ij[1]*r1_kj[0]; /* 9 */
1700
1701 b1[0] = r1_kl[1]*r1_kj[2]-r1_kl[2]*r1_kj[1];
1702 b1[1] = r1_kl[2]*r1_kj[0]-r1_kl[0]*r1_kj[2];
1703 b1[2] = r1_kl[0]*r1_kj[1]-r1_kl[1]*r1_kj[0]; /* 9 */
1704
1705 tmp = pbc_rvec_sub(pbc,x[a1l],x[a1k],h1);
1706
1707 ra21 = iprod(a1,a1); /* 5 */
1708 rb21 = iprod(b1,b1); /* 5 */
1709 rg21 = iprod(r1_kj,r1_kj); /* 5 */
1710 rg1 = sqrt(rg21);
1711
1712 rgr1 = 1.0/rg1;
1713 ra2r1 = 1.0/ra21;
1714 rb2r1 = 1.0/rb21;
1715 rabr1 = sqrt(ra2r1*rb2r1);
1716
1717 sin_phi1 = rg1 * rabr1 * iprod(a1,h1) * (-1);
1718
1719 if(cos_phi1 < -0.5 || cos_phi1 > 0.5)
1720 {
1721 phi1 = asin(sin_phi1);
1722
1723 if(cos_phi1 < 0)
1724 {
1725 if(phi1 > 0)
1726 {
1727 phi1 = M_PI - phi1;
1728 }
1729 else
1730 {
1731 phi1 = -M_PI - phi1;
1732 }
1733 }
1734 }
1735 else
1736 {
1737 phi1 = acos(cos_phi1);
1738
1739 if(sin_phi1 < 0)
1740 {
1741 phi1 = -phi1;
1742 }
1743 }
1744
1745 xphi1 = phi1 + M_PI; /* 1 */
1746
1747 /* Second torsion */
1748 a2i = aj;
1749 a2j = ak;
1750 a2k = al;
1751 a2l = am;
1752
1753 phi2 = dih_angle(x[a2i], x[a2j], x[a2k], x[a2l], pbc, r2_ij, r2_kj, r2_kl, m2, n2,
1754 &sign2, &t12, &t22, &t32); /* 84 */
1755
1756 cos_phi2 = cos(phi2);
1757
1758 a2[0] = r2_ij[1]*r2_kj[2]-r2_ij[2]*r2_kj[1];
1759 a2[1] = r2_ij[2]*r2_kj[0]-r2_ij[0]*r2_kj[2];
1760 a2[2] = r2_ij[0]*r2_kj[1]-r2_ij[1]*r2_kj[0]; /* 9 */
1761
1762 b2[0] = r2_kl[1]*r2_kj[2]-r2_kl[2]*r2_kj[1];
1763 b2[1] = r2_kl[2]*r2_kj[0]-r2_kl[0]*r2_kj[2];
1764 b2[2] = r2_kl[0]*r2_kj[1]-r2_kl[1]*r2_kj[0]; /* 9 */
1765
1766 tmp = pbc_rvec_sub(pbc,x[a2l],x[a2k],h2);
1767
1768 ra22 = iprod(a2,a2); /* 5 */
1769 rb22 = iprod(b2,b2); /* 5 */
1770 rg22 = iprod(r2_kj,r2_kj); /* 5 */
1771 rg2 = sqrt(rg22);
1772
1773 rgr2 = 1.0/rg2;
1774 ra2r2 = 1.0/ra22;
1775 rb2r2 = 1.0/rb22;
1776 rabr2 = sqrt(ra2r2*rb2r2);
1777
1778 sin_phi2 = rg2 * rabr2 * iprod(a2,h2) * (-1);
1779
1780 if(cos_phi2 < -0.5 || cos_phi2 > 0.5)
1781 {
1782 phi2 = asin(sin_phi2);
1783
1784 if(cos_phi2 < 0)
1785 {
1786 if(phi2 > 0)
1787 {
1788 phi2 = M_PI - phi2;
1789 }
1790 else
1791 {
1792 phi2 = -M_PI - phi2;
1793 }
1794 }
1795 }
1796 else
1797 {
1798 phi2 = acos(cos_phi2);
1799
1800 if(sin_phi2 < 0)
1801 {
1802 phi2 = -phi2;
1803 }
1804 }
1805
1806 xphi2 = phi2 + M_PI; /* 1 */
1807
1808 /* Range mangling */
1809 if(xphi1<0)
1810 {
1811 xphi1 = xphi1 + 2*M_PI;
1812 }
1813 else if(xphi1>=2*M_PI)
1814 {
1815 xphi1 = xphi1 - 2*M_PI;
1816 }
1817
1818 if(xphi2<0)
1819 {
1820 xphi2 = xphi2 + 2*M_PI;
1821 }
1822 else if(xphi2>=2*M_PI)
1823 {
1824 xphi2 = xphi2 - 2*M_PI;
1825 }
1826
1827 /* Number of grid points */
1828 dx = 2*M_PI / cmap_grid->grid_spacing;
1829
1830 /* Where on the grid are we */
1831 iphi1 = (int)(xphi1/dx);
1832 iphi2 = (int)(xphi2/dx);
1833
1834 iphi1 = cmap_setup_grid_index(iphi1, cmap_grid->grid_spacing, &ip1m1,&ip1p1,&ip1p2);
1835 iphi2 = cmap_setup_grid_index(iphi2, cmap_grid->grid_spacing, &ip2m1,&ip2p1,&ip2p2);
1836
1837 pos1 = iphi1*cmap_grid->grid_spacing+iphi2;
1838 pos2 = ip1p1*cmap_grid->grid_spacing+iphi2;
1839 pos3 = ip1p1*cmap_grid->grid_spacing+ip2p1;
1840 pos4 = iphi1*cmap_grid->grid_spacing+ip2p1;
1841
1842 ty[0] = cmapd[pos1*4];
1843 ty[1] = cmapd[pos2*4];
1844 ty[2] = cmapd[pos3*4];
1845 ty[3] = cmapd[pos4*4];
1846
1847 ty1[0] = cmapd[pos1*4+1];
1848 ty1[1] = cmapd[pos2*4+1];
1849 ty1[2] = cmapd[pos3*4+1];
1850 ty1[3] = cmapd[pos4*4+1];
1851
1852 ty2[0] = cmapd[pos1*4+2];
1853 ty2[1] = cmapd[pos2*4+2];
1854 ty2[2] = cmapd[pos3*4+2];
1855 ty2[3] = cmapd[pos4*4+2];
1856
1857 ty12[0] = cmapd[pos1*4+3];
1858 ty12[1] = cmapd[pos2*4+3];
1859 ty12[2] = cmapd[pos3*4+3];
1860 ty12[3] = cmapd[pos4*4+3];
1861
1862 /* Switch to degrees */
1863 dx = 360.0 / cmap_grid->grid_spacing;
1864 xphi1 = xphi1 * RAD2DEG;
1865 xphi2 = xphi2 * RAD2DEG;
1866
1867 for(i=0;i<4;i++) /* 16 */
1868 {
1869 tx[i] = ty[i];
1870 tx[i+4] = ty1[i]*dx;
1871 tx[i+8] = ty2[i]*dx;
1872 tx[i+12] = ty12[i]*dx*dx;
1873 }
1874
1875 idx=0;
1876 for(i=0;i<4;i++) /* 1056 */
1877 {
1878 for(j=0;j<4;j++)
1879 {
1880 xx = 0;
1881 for(k=0;k<16;k++)
1882 {
1883 xx = xx + cmap_coeff_matrix[k*16+idx]*tx[k];
1884 }
1885
1886 idx++;
1887 tc[i*4+j]=xx;
1888 }
1889 }
1890
1891 tt = (xphi1-iphi1*dx)/dx;
1892 tu = (xphi2-iphi2*dx)/dx;
1893
1894 e = 0;
1895 df1 = 0;
1896 df2 = 0;
1897 ddf1 = 0;
1898 ddf2 = 0;
1899 ddf12 = 0;
1900
1901 for(i=3;i>=0;i--)
1902 {
1903 l1 = loop_index[i][3];
1904 l2 = loop_index[i][2];
1905 l3 = loop_index[i][1];
1906
1907 e = tt * e + ((tc[i*4+3]*tu+tc[i*4+2])*tu + tc[i*4+1])*tu+tc[i*4];
1908 df1 = tu * df1 + (3.0*tc[l1]*tt+2.0*tc[l2])*tt+tc[l3];
1909 df2 = tt * df2 + (3.0*tc[i*4+3]*tu+2.0*tc[i*4+2])*tu+tc[i*4+1];
1910 ddf1 = tu * ddf1 + 2.0*3.0*tc[l1]*tt+2.0*tc[l2];
1911 ddf2 = tt * ddf2 + 2.0*3.0*tc[4*i+3]*tu+2.0*tc[4*i+2];
1912 }
1913
1914 ddf12 = tc[5] + 2.0*tc[9]*tt + 3.0*tc[13]*tt*tt + 2.0*tu*(tc[6]+2.0*tc[10]*tt+3.0*tc[14]*tt*tt) +
1915 3.0*tu*tu*(tc[7]+2.0*tc[11]*tt+3.0*tc[15]*tt*tt);
1916
1917 fac = RAD2DEG/dx;
1918 df1 = df1 * fac;
1919 df2 = df2 * fac;
1920 ddf1 = ddf1 * fac * fac;
1921 ddf2 = ddf2 * fac * fac;
1922 ddf12 = ddf12 * fac * fac;
1923
1924 /* CMAP energy */
1925 vtot += e;
1926
1927 /* Do forces - first torsion */
1928 fg1 = iprod(r1_ij,r1_kj);
1929 hg1 = iprod(r1_kl,r1_kj);
1930 fga1 = fg1*ra2r1*rgr1;
1931 hgb1 = hg1*rb2r1*rgr1;
1932 gaa1 = -ra2r1*rg1;
1933 gbb1 = rb2r1*rg1;
1934
1935 for(i=0;i<DIM;i++)
1936 {
1937 dtf1[i] = gaa1 * a1[i];
1938 dtg1[i] = fga1 * a1[i] - hgb1 * b1[i];
1939 dth1[i] = gbb1 * b1[i];
1940
1941 f1[i] = df1 * dtf1[i];
1942 g1[i] = df1 * dtg1[i];
1943 h1[i] = df1 * dth1[i];
1944
1945 f1_i[i] = f1[i];
1946 f1_j[i] = -f1[i] - g1[i];
1947 f1_k[i] = h1[i] + g1[i];
1948 f1_l[i] = -h1[i];
1949
1950 f[a1i][i] = f[a1i][i] + f1_i[i];
1951 f[a1j][i] = f[a1j][i] + f1_j[i]; /* - f1[i] - g1[i] */
1952 f[a1k][i] = f[a1k][i] + f1_k[i]; /* h1[i] + g1[i] */
1953 f[a1l][i] = f[a1l][i] + f1_l[i]; /* h1[i] */
1954 }
1955
1956 /* Do forces - second torsion */
1957 fg2 = iprod(r2_ij,r2_kj);
1958 hg2 = iprod(r2_kl,r2_kj);
1959 fga2 = fg2*ra2r2*rgr2;
1960 hgb2 = hg2*rb2r2*rgr2;
1961 gaa2 = -ra2r2*rg2;
1962 gbb2 = rb2r2*rg2;
1963
1964 for(i=0;i<DIM;i++)
1965 {
1966 dtf2[i] = gaa2 * a2[i];
1967 dtg2[i] = fga2 * a2[i] - hgb2 * b2[i];
1968 dth2[i] = gbb2 * b2[i];
1969
1970 f2[i] = df2 * dtf2[i];
1971 g2[i] = df2 * dtg2[i];
1972 h2[i] = df2 * dth2[i];
1973
1974 f2_i[i] = f2[i];
1975 f2_j[i] = -f2[i] - g2[i];
1976 f2_k[i] = h2[i] + g2[i];
1977 f2_l[i] = -h2[i];
1978
1979 f[a2i][i] = f[a2i][i] + f2_i[i]; /* f2[i] */
1980 f[a2j][i] = f[a2j][i] + f2_j[i]; /* - f2[i] - g2[i] */
1981 f[a2k][i] = f[a2k][i] + f2_k[i]; /* h2[i] + g2[i] */
1982 f[a2l][i] = f[a2l][i] + f2_l[i]; /* - h2[i] */
1983 }
1984
1985 /* Shift forces */
1986 if(g)
1987 {
1988 copy_ivec(SHIFT_IVEC(g,a1j), jt1);
1989 ivec_sub(SHIFT_IVEC(g,a1i), jt1,dt1_ij);
1990 ivec_sub(SHIFT_IVEC(g,a1k), jt1,dt1_kj);
1991 ivec_sub(SHIFT_IVEC(g,a1l), jt1,dt1_lj);
1992 t11 = IVEC2IS(dt1_ij);
1993 t21 = IVEC2IS(dt1_kj);
1994 t31 = IVEC2IS(dt1_lj);
1995
1996 copy_ivec(SHIFT_IVEC(g,a2j), jt2);
1997 ivec_sub(SHIFT_IVEC(g,a2i), jt2,dt2_ij);
1998 ivec_sub(SHIFT_IVEC(g,a2k), jt2,dt2_kj);
1999 ivec_sub(SHIFT_IVEC(g,a2l), jt2,dt2_lj);
2000 t12 = IVEC2IS(dt2_ij);
2001 t22 = IVEC2IS(dt2_kj);
2002 t32 = IVEC2IS(dt2_lj);
2003 }
2004 else if(pbc)
2005 {
2006 t31 = pbc_rvec_sub(pbc,x[a1l],x[a1j],h1);
2007 t32 = pbc_rvec_sub(pbc,x[a2l],x[a2j],h2);
2008 }
2009 else
2010 {
2011 t31 = CENTRAL;
2012 t32 = CENTRAL;
2013 }
2014
2015 rvec_inc(fshift[t11],f1_i);
2016 rvec_inc(fshift[CENTRAL],f1_j);
2017 rvec_inc(fshift[t21],f1_k);
2018 rvec_inc(fshift[t31],f1_l);
2019
2020 rvec_inc(fshift[t21],f2_i);
2021 rvec_inc(fshift[CENTRAL],f2_j);
2022 rvec_inc(fshift[t22],f2_k);
2023 rvec_inc(fshift[t32],f2_l);
2024 }
2025 return vtot;
2026 }
2027
2028
2029
2030 /***********************************************************
2031 *
2032 * G R O M O S 9 6 F U N C T I O N S
2033 *
2034 ***********************************************************/
2035 real g96harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
2036 real *V,real *F)
2037 {
2038 const real half=0.5;
2039 real L1,kk,x0,dx,dx2;
2040 real v,f,dvdl;
2041
2042 L1 = 1.0-lambda;
2043 kk = L1*kA+lambda*kB;
2044 x0 = L1*xA+lambda*xB;
2045
2046 dx = x-x0;
2047 dx2 = dx*dx;
2048
2049 f = -kk*dx;
2050 v = half*kk*dx2;
2051 dvdl = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
2052
2053 *F = f;
2054 *V = v;
2055
2056 return dvdl;
2057
2058 /* That was 21 flops */
2059 }
2060
2061 real g96bonds(int nbonds,
2062 const t_iatom forceatoms[],const t_iparams forceparams[],
2063 const rvec x[],rvec f[],rvec fshift[],
2064 const t_pbc *pbc,const t_graph *g,
2065 real lambda,real *dvdlambda,
2066 const t_mdatoms *md,t_fcdata *fcd,
2067 int *global_atom_index)
2068 {
2069 int i,m,ki,ai,aj,type;
2070 real dr2,fbond,vbond,fij,vtot;
2071 rvec dx;
2072 ivec dt;
2073
2074 vtot = 0.0;
2075 for(i=0; (i<nbonds); ) {
2076 type = forceatoms[i++];
2077 ai = forceatoms[i++];
2078 aj = forceatoms[i++];
2079
2080 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
2081 dr2 = iprod(dx,dx); /* 5 */
2082
2083 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2084 forceparams[type].harmonic.krB,
2085 forceparams[type].harmonic.rA,
2086 forceparams[type].harmonic.rB,
2087 dr2,lambda,&vbond,&fbond);
2088
2089 vtot += 0.5*vbond; /* 1*/
2090 #ifdef DEBUG
2091 if (debug)
2092 fprintf(debug,"G96-BONDS: dr = %10g vbond = %10g fbond = %10g\n",
2093 sqrt(dr2),vbond,fbond);
2094 #endif
2095
2096 if (g) {
2097 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2098 ki=IVEC2IS(dt);
2099 }
2100 for (m=0; (m<DIM); m++) { /* 15 */
2101 fij=fbond*dx[m];
2102 f[ai][m]+=fij;
2103 f[aj][m]-=fij;
2104 fshift[ki][m]+=fij;
2105 fshift[CENTRAL][m]-=fij;
2106 }
2107 } /* 44 TOTAL */
2108 return vtot;
2109 }
2110
2111 real g96bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
2112 rvec r_ij,rvec r_kj,
2113 int *t1,int *t2)
2114 /* Return value is the angle between the bonds i-j and j-k */
2115 {
2116 real costh;
2117
2118 *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij); /* 3 */
2119 *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj); /* 3 */
2120
2121 costh=cos_angle(r_ij,r_kj); /* 25 */
2122 /* 41 TOTAL */
2123 return costh;
2124 }
2125
2126 real g96angles(int nbonds,
2127 const t_iatom forceatoms[],const t_iparams forceparams[],
2128 const rvec x[],rvec f[],rvec fshift[],
2129 const t_pbc *pbc,const t_graph *g,
2130 real lambda,real *dvdlambda,
2131 const t_mdatoms *md,t_fcdata *fcd,
2132 int *global_atom_index)
2133 {
2134 int i,ai,aj,ak,type,m,t1,t2;
2135 rvec r_ij,r_kj;
2136 real cos_theta,dVdt,va,vtot;
2137 real rij_1,rij_2,rkj_1,rkj_2,rijrkj_1;
2138 rvec f_i,f_j,f_k;
2139 ivec jt,dt_ij,dt_kj;
2140
2141 vtot = 0.0;
2142 for(i=0; (i<nbonds); ) {
2143 type = forceatoms[i++];
2144 ai = forceatoms[i++];
2145 aj = forceatoms[i++];
2146 ak = forceatoms[i++];
2147
2148 cos_theta = g96bond_angle(x[ai],x[aj],x[ak],pbc,r_ij,r_kj,&t1,&t2);
2149
2150 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
2151 forceparams[type].harmonic.krB,
2152 forceparams[type].harmonic.rA,
2153 forceparams[type].harmonic.rB,
2154 cos_theta,lambda,&va,&dVdt);
2155 vtot += va;
2156
2157 rij_1 = gmx_invsqrt(iprod(r_ij,r_ij));
2158 rkj_1 = gmx_invsqrt(iprod(r_kj,r_kj));
2159 rij_2 = rij_1*rij_1;
2160 rkj_2 = rkj_1*rkj_1;
2161 rijrkj_1 = rij_1*rkj_1; /* 23 */
2162
2163 #ifdef DEBUG
2164 if (debug)
2165 fprintf(debug,"G96ANGLES: costheta = %10g vth = %10g dV/dct = %10g\n",
2166 cos_theta,va,dVdt);
2167 #endif
2168 for (m=0; (m<DIM); m++) { /* 42 */
2169 f_i[m]=dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
2170 f_k[m]=dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
2171 f_j[m]=-f_i[m]-f_k[m];
2172 f[ai][m]+=f_i[m];
2173 f[aj][m]+=f_j[m];
2174 f[ak][m]+=f_k[m];
2175 }
2176
2177 if (g) {
2178 copy_ivec(SHIFT_IVEC(g,aj),jt);
2179
2180 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2181 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2182 t1=IVEC2IS(dt_ij);
2183 t2=IVEC2IS(dt_kj);
2184 }
2185 rvec_inc(fshift[t1],f_i);
2186 rvec_inc(fshift[CENTRAL],f_j);
2187 rvec_inc(fshift[t2],f_k); /* 9 */
2188 /* 163 TOTAL */
2189 }
2190 return vtot;
2191 }
2192
2193 real cross_bond_bond(int nbonds,
2194 const t_iatom forceatoms[],const t_iparams forceparams[],
2195 const rvec x[],rvec f[],rvec fshift[],
2196 const t_pbc *pbc,const t_graph *g,
2197 real lambda,real *dvdlambda,
2198 const t_mdatoms *md,t_fcdata *fcd,
2199 int *global_atom_index)
2200 {
2201 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2202 * pp. 842-847
2203 */
2204 int i,ai,aj,ak,type,m,t1,t2;
2205 rvec r_ij,r_kj;
2206 real vtot,vrr,s1,s2,r1,r2,r1e,r2e,krr;
2207 rvec f_i,f_j,f_k;
2208 ivec jt,dt_ij,dt_kj;
2209
2210 vtot = 0.0;
2211 for(i=0; (i<nbonds); ) {
2212 type = forceatoms[i++];
2213 ai = forceatoms[i++];
2214 aj = forceatoms[i++];
2215 ak = forceatoms[i++];
2216 r1e = forceparams[type].cross_bb.r1e;
2217 r2e = forceparams[type].cross_bb.r2e;
2218 krr = forceparams[type].cross_bb.krr;
2219
2220 /* Compute distance vectors ... */
2221 t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2222 t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2223
2224 /* ... and their lengths */
2225 r1 = norm(r_ij);
2226 r2 = norm(r_kj);
2227
2228 /* Deviations from ideality */
2229 s1 = r1-r1e;
2230 s2 = r2-r2e;
2231
2232 /* Energy (can be negative!) */
2233 vrr = krr*s1*s2;
2234 vtot += vrr;
2235
2236 /* Forces */
2237 svmul(-krr*s2/r1,r_ij,f_i);
2238 svmul(-krr*s1/r2,r_kj,f_k);
2239
2240 for (m=0; (m<DIM); m++) { /* 12 */
2241 f_j[m] = -f_i[m] - f_k[m];
2242 f[ai][m] += f_i[m];
2243 f[aj][m] += f_j[m];
2244 f[ak][m] += f_k[m];
2245 }
2246
2247 /* Virial stuff */
2248 if (g) {
2249 copy_ivec(SHIFT_IVEC(g,aj),jt);
2250
2251 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2252 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2253 t1=IVEC2IS(dt_ij);
2254 t2=IVEC2IS(dt_kj);
2255 }
2256 rvec_inc(fshift[t1],f_i);
2257 rvec_inc(fshift[CENTRAL],f_j);
2258 rvec_inc(fshift[t2],f_k); /* 9 */
2259 /* 163 TOTAL */
2260 }
2261 return vtot;
2262 }
2263
2264 real cross_bond_angle(int nbonds,
2265 const t_iatom forceatoms[],const t_iparams forceparams[],
2266 const rvec x[],rvec f[],rvec fshift[],
2267 const t_pbc *pbc,const t_graph *g,
2268 real lambda,real *dvdlambda,
2269 const t_mdatoms *md,t_fcdata *fcd,
2270 int *global_atom_index)
2271 {
2272 /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
2273 * pp. 842-847
2274 */
2275 int i,ai,aj,ak,type,m,t1,t2,t3;
2276 rvec r_ij,r_kj,r_ik;
2277 real vtot,vrt,s1,s2,s3,r1,r2,r3,r1e,r2e,r3e,krt,k1,k2,k3;
2278 rvec f_i,f_j,f_k;
2279 ivec jt,dt_ij,dt_kj;
2280
2281 vtot = 0.0;
2282 for(i=0; (i<nbonds); ) {
2283 type = forceatoms[i++];
2284 ai = forceatoms[i++];
2285 aj = forceatoms[i++];
2286 ak = forceatoms[i++];
2287 r1e = forceparams[type].cross_ba.r1e;
2288 r2e = forceparams[type].cross_ba.r2e;
2289 r3e = forceparams[type].cross_ba.r3e;
2290 krt = forceparams[type].cross_ba.krt;
2291
2292 /* Compute distance vectors ... */
2293 t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
2294 t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
2295 t3 = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);
2296
2297 /* ... and their lengths */
2298 r1 = norm(r_ij);
2299 r2 = norm(r_kj);
2300 r3 = norm(r_ik);
2301
2302 /* Deviations from ideality */
2303 s1 = r1-r1e;
2304 s2 = r2-r2e;
2305 s3 = r3-r3e;
2306
2307 /* Energy (can be negative!) */
2308 vrt = krt*s3*(s1+s2);
2309 vtot += vrt;
2310
2311 /* Forces */
2312 k1 = -krt*(s3/r1);
2313 k2 = -krt*(s3/r2);
2314 k3 = -krt*(s1+s2)/r3;
2315 for(m=0; (m<DIM); m++) {
2316 f_i[m] = k1*r_ij[m] + k3*r_ik[m];
2317 f_k[m] = k2*r_kj[m] - k3*r_ik[m];
2318 f_j[m] = -f_i[m] - f_k[m];
2319 }
2320
2321 for (m=0; (m<DIM); m++) { /* 12 */
2322 f[ai][m] += f_i[m];
2323 f[aj][m] += f_j[m];
2324 f[ak][m] += f_k[m];
2325 }
2326
2327 /* Virial stuff */
2328 if (g) {
2329 copy_ivec(SHIFT_IVEC(g,aj),jt);
2330
2331 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2332 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2333 t1=IVEC2IS(dt_ij);
2334 t2=IVEC2IS(dt_kj);
2335 }
2336 rvec_inc(fshift[t1],f_i);
2337 rvec_inc(fshift[CENTRAL],f_j);
2338 rvec_inc(fshift[t2],f_k); /* 9 */
2339 /* 163 TOTAL */
2340 }
2341 return vtot;
2342 }
2343
2344 static real bonded_tab(const char *type,int table_nr,
2345 const bondedtable_t *table,real kA,real kB,real r,
2346 real lambda,real *V,real *F)
2347 {
2348 real k,tabscale,*VFtab,rt,eps,eps2,Yt,Ft,Geps,Heps2,Fp,VV,FF;
2349 int n0,nnn;
2350 real v,f,dvdl;
2351
2352 k = (1.0 - lambda)*kA + lambda*kB;
2353
2354 tabscale = table->scale;
2355 VFtab = table->tab;
2356
2357 rt = r*tabscale;
2358 n0 = rt;
2359 if (n0 >= table->n) {
2360 gmx_fatal(FARGS,"A tabulated %s interaction table number %d is out of the table range: r %f, between table indices %d and %d, table length %d",
2361 type,table_nr,r,n0,n0+1,table->n);
2362 }
2363 eps = rt - n0;
2364 eps2 = eps*eps;
2365 nnn = 4*n0;
2366 Yt = VFtab[nnn];
2367 Ft = VFtab[nnn+1];
2368 Geps = VFtab[nnn+2]*eps;
2369 Heps2 = VFtab[nnn+3]*eps2;
2370 Fp = Ft + Geps + Heps2;
2371 VV = Yt + Fp*eps;
2372 FF = Fp + Geps + 2.0*Heps2;
2373
2374 *F = -k*FF*tabscale;
2375 *V = k*VV;
2376 dvdl = (kB - kA)*VV;
2377
2378 return dvdl;
2379
2380 /* That was 22 flops */
2381 }
2382
2383 real tab_bonds(int nbonds,
2384 const t_iatom forceatoms[],const t_iparams forceparams[],
2385 const rvec x[],rvec f[],rvec fshift[],
2386 const t_pbc *pbc,const t_graph *g,
2387 real lambda,real *dvdlambda,
2388 const t_mdatoms *md,t_fcdata *fcd,
2389 int *global_atom_index)
2390 {
2391 int i,m,ki,ai,aj,type,table;
2392 real dr,dr2,fbond,vbond,fij,vtot;
2393 rvec dx;
2394 ivec dt;
2395
2396 vtot = 0.0;
2397 for(i=0; (i<nbonds); ) {
2398 type = forceatoms[i++];
2399 ai = forceatoms[i++];
2400 aj = forceatoms[i++];
2401
2402 ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx); /* 3 */
2403 dr2 = iprod(dx,dx); /* 5 */
2404 dr = dr2*gmx_invsqrt(dr2); /* 10 */
2405
2406 table = forceparams[type].tab.table;
2407
2408 *dvdlambda += bonded_tab("bond",table,
2409 &fcd->bondtab[table],
2410 forceparams[type].tab.kA,
2411 forceparams[type].tab.kB,
2412 dr,lambda,&vbond,&fbond); /* 22 */
2413
2414 if (dr2 == 0.0)
2415 continue;
2416
2417
2418 vtot += vbond;/* 1*/
2419 fbond *= gmx_invsqrt(dr2); /* 6 */
2420 #ifdef DEBUG
2421 if (debug)
2422 fprintf(debug,"TABBONDS: dr = %10g vbond = %10g fbond = %10g\n",
2423 dr,vbond,fbond);
2424 #endif
2425 if (g) {
2426 ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
2427 ki=IVEC2IS(dt);
2428 }
2429 for (m=0; (m<DIM); m++) { /* 15 */
2430 fij=fbond*dx[m];
2431 f[ai][m]+=fij;
2432 f[aj][m]-=fij;
2433 fshift[ki][m]+=fij;
2434 fshift[CENTRAL][m]-=fij;
2435 }
2436 } /* 62 TOTAL */
2437 return vtot;
2438 }
2439
2440 real tab_angles(int nbonds,
2441 const t_iatom forceatoms[],const t_iparams forceparams[],
2442 const rvec x[],rvec f[],rvec fshift[],
2443 const t_pbc *pbc,const t_graph *g,
2444 real lambda,real *dvdlambda,
2445 const t_mdatoms *md,t_fcdata *fcd,
2446 int *global_atom_index)
2447 {
2448 int i,ai,aj,ak,t1,t2,type,table;
2449 rvec r_ij,r_kj;
2450 real cos_theta,cos_theta2,theta,dVdt,va,vtot;
2451 ivec jt,dt_ij,dt_kj;
2452
2453 vtot = 0.0;
2454 for(i=0; (i<nbonds); ) {
2455 type = forceatoms[i++];
2456 ai = forceatoms[i++];
2457 aj = forceatoms[i++];
2458 ak = forceatoms[i++];
2459
2460 theta = bond_angle(x[ai],x[aj],x[ak],pbc,
2461 r_ij,r_kj,&cos_theta,&t1,&t2); /* 41 */
2462
2463 table = forceparams[type].tab.table;
2464
2465 *dvdlambda += bonded_tab("angle",table,
2466 &fcd->angletab[table],
2467 forceparams[type].tab.kA,
2468 forceparams[type].tab.kB,
2469 theta,lambda,&va,&dVdt); /* 22 */
2470 vtot += va;
2471
2472 cos_theta2 = sqr(cos_theta); /* 1 */
2473 if (cos_theta2 < 1) {
2474 int m;
2475 real snt,st,sth;
2476 real cik,cii,ckk;
2477 real nrkj2,nrij2;
2478 rvec f_i,f_j,f_k;
2479
2480 st = dVdt*gmx_invsqrt(1 - cos_theta2); /* 12 */
2481 sth = st*cos_theta; /* 1 */
2482 #ifdef DEBUG
2483 if (debug)
2484 fprintf(debug,"ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n",
2485 theta*RAD2DEG,va,dVdt);
2486 #endif
2487 nrkj2=iprod(r_kj,r_kj); /* 5 */
2488 nrij2=iprod(r_ij,r_ij);
2489
2490 cik=st*gmx_invsqrt(nrkj2*nrij2); /* 12 */
2491 cii=sth/nrij2; /* 10 */
2492 ckk=sth/nrkj2; /* 10 */
2493
2494 for (m=0; (m<DIM); m++) { /* 39 */
2495 f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
2496 f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
2497 f_j[m]=-f_i[m]-f_k[m];
2498 f[ai][m]+=f_i[m];
2499 f[aj][m]+=f_j[m];
2500 f[ak][m]+=f_k[m];
2501 }
2502 if (g) {
2503 copy_ivec(SHIFT_IVEC(g,aj),jt);
2504
2505 ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
2506 ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
2507 t1=IVEC2IS(dt_ij);
2508 t2=IVEC2IS(dt_kj);
2509 }
2510 rvec_inc(fshift[t1],f_i);
2511 rvec_inc(fshift[CENTRAL],f_j);
2512 rvec_inc(fshift[t2],f_k);
2513 } /* 169 TOTAL */
2514 }
2515 return vtot;
2516 }
2517
2518 real tab_dihs(int nbonds,
2519 const t_iatom forceatoms[],const t_iparams forceparams[],
2520 const rvec x[],rvec f[],rvec fshift[],
2521 const t_pbc *pbc,const t_graph *g,
2522 real lambda,real *dvdlambda,
2523 const t_mdatoms *md,t_fcdata *fcd,
2524 int *global_atom_index)
2525 {
2526 int i,type,ai,aj,ak,al,table;
2527 int t1,t2,t3;
2528 rvec r_ij,r_kj,r_kl,m,n;
2529 real phi,sign,ddphi,vpd,vtot;
2530
2531 vtot = 0.0;
2532 for(i=0; (i<nbonds); ) {
2533 type = forceatoms[i++];
2534 ai = forceatoms[i++];
2535 aj = forceatoms[i++];
2536 ak = forceatoms[i++];
2537 al = forceatoms[i++];
2538
2539 phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
2540 &sign,&t1,&t2,&t3); /* 84 */
2541
2542 table = forceparams[type].tab.table;
2543
2544 /* Hopefully phi+M_PI never results in values < 0 */
2545 *dvdlambda += bonded_tab("dihedral",table,
2546 &fcd->dihtab[table],
2547 forceparams[type].tab.kA,
2548 forceparams[type].tab.kB,
2549 phi+M_PI,lambda,&vpd,&ddphi);
2550
2551 vtot += vpd;
2552 do_dih_fup(ai,aj,ak,al,-ddphi,r_ij,r_kj,r_kl,m,n,
2553 f,fshift,pbc,g,x,t1,t2,t3); /* 112 */
2554
2555 #ifdef DEBUG
2556 fprintf(debug,"pdih: (%d,%d,%d,%d) phi=%g\n",
2557 ai,aj,ak,al,phi);
2558 #endif
2559 } /* 227 TOTAL */
2560
2561 return vtot;
2562 }
2563
2564 void calc_bonds(FILE *fplog,const gmx_multisim_t *ms,
2565 const t_idef *idef,
2566 rvec x[],history_t *hist,
2567 rvec f[],t_forcerec *fr,
2568 const t_pbc *pbc,const t_graph *g,
2569 gmx_enerdata_t *enerd,t_nrnb *nrnb,
2570 real lambda,
2571 const t_mdatoms *md,
2572 t_fcdata *fcd,int *global_atom_index,
2573 t_atomtypes *atype, gmx_genborn_t *born,
2574 gmx_bool bPrintSepPot,gmx_large_int_t step)
2575 {
2576 int ftype,nbonds,ind,nat1;
2577 real *epot,v,dvdl;
2578 const t_pbc *pbc_null;
2579 char buf[22];
2580
2581 if (fr->bMolPBC)
2582 pbc_null = pbc;
2583 else
2584 pbc_null = NULL;
2585
2586 if (bPrintSepPot)
2587 fprintf(fplog,"Step %s: bonded V and dVdl for this node\n",
2588 gmx_step_str(step,buf));
2589
2590 #ifdef DEBUG
2591 if (g && debug)
2592 p_graph(debug,"Bondage is fun",g);
2593 #endif
2594
2595 epot = enerd->term;
2596
2597 /* Do pre force calculation stuff which might require communication */
2598 if (idef->il[F_ORIRES].nr) {
2599 epot[F_ORIRESDEV] = calc_orires_dev(ms,idef->il[F_ORIRES].nr,
2600 idef->il[F_ORIRES].iatoms,
2601 idef->iparams,md,(const rvec*)x,
2602 pbc_null,fcd,hist);
2603 }
2604 if (idef->il[F_DISRES].nr) {
2605 calc_disres_R_6(ms,idef->il[F_DISRES].nr,
2606 idef->il[F_DISRES].iatoms,
2607 idef->iparams,(const rvec*)x,pbc_null,
2608 fcd,hist);
2609 }
2610
2611 /* Loop over all bonded force types to calculate the bonded forces */
2612 for(ftype=0; (ftype<F_NRE); ftype++) {
2613 if(ftype<F_GB12 || ftype>F_GB14) {
2614 if (interaction_function[ftype].flags & IF_BOND &&
2615 !(ftype == F_CONNBONDS || ftype == F_POSRES)) {
2616 nbonds=idef->il[ftype].nr;
2617 if (nbonds > 0) {
2618 ind = interaction_function[ftype].nrnb_ind;
2619 nat1 = interaction_function[ftype].nratoms + 1;
2620 dvdl = 0;
2621 if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB) {
2622 if(ftype==F_CMAP)
2623 {
2624 v = cmap_dihs(nbonds,idef->il[ftype].iatoms,
2625 idef->iparams,&idef->cmap_grid,
2626 (const rvec*)x,f,fr->fshift,
2627 pbc_null,g,lambda,&dvdl,md,fcd,
2628 global_atom_index);
2629 }
2630 else
2631 {
2632 v =
2633 interaction_function[ftype].ifunc(nbonds,idef->il[ftype].iatoms,
2634 idef->iparams,
2635 (const rvec*)x,f,fr->fshift,
2636 pbc_null,g,lambda,&dvdl,md,fcd,
2637 global_atom_index);
2638 }
2639
2640 if (bPrintSepPot) {
2641 fprintf(fplog," %-23s #%4d V %12.5e dVdl %12.5e\n",
2642 interaction_function[ftype].longname,nbonds/nat1,v,dvdl);
2643 }
2644 } else {
2645 v = do_listed_vdw_q(ftype,nbonds,idef->il[ftype].iatoms,
2646 idef->iparams,
2647 (const rvec*)x,f,fr->fshift,
2648 pbc_null,g,
2649 lambda,&dvdl,
2650 md,fr,&enerd->grpp,global_atom_index);
2651 if (bPrintSepPot) {
2652 fprintf(fplog," %-5s + %-15s #%4d dVdl %12.5e\n",
2653 interaction_function[ftype].longname,
2654 interaction_function[F_COUL14].longname,nbonds/nat1,dvdl);
2655 }
2656 }
2657 if (ind != -1)
2658 inc_nrnb(nrnb,ind,nbonds/nat1);
2659 epot[ftype] += v;
2660 enerd->dvdl_nonlin += dvdl;
2661 }
2662 }
2663 }
2664 }
2665 /* Copy the sum of violations for the distance restraints from fcd */
2666 if (fcd)
2667 epot[F_DISRESVIOL] = fcd->disres.sumviol;
2668 }
2669
2670 void calc_bonds_lambda(FILE *fplog,
2671 const t_idef *idef,
2672 rvec x[],
2673 t_forcerec *fr,
2674 const t_pbc *pbc,const t_graph *g,
2675 gmx_enerdata_t *enerd,t_nrnb *nrnb,
2676 real lambda,
2677 const t_mdatoms *md,
2678 t_fcdata *fcd,int *global_atom_index)
2679 {
2680 int ftype,nbonds_np,nbonds,ind, nat1;
2681 real *epot,v,dvdl;
2682 rvec *f,*fshift_orig;
2683 const t_pbc *pbc_null;
2684 t_iatom *iatom_fe;
2685
2686 if (fr->bMolPBC)
2687 pbc_null = pbc;
2688 else
2689 pbc_null = NULL;
2690
2691 epot = enerd->term;
2692
2693 snew(f,fr->natoms_force);
2694 /* We want to preserve the fshift array in forcerec */
2695 fshift_orig = fr->fshift;
2696 snew(fr->fshift,SHIFTS);
2697
2698 /* Loop over all bonded force types to calculate the bonded forces */
2699 for(ftype=0; (ftype<F_NRE); ftype++) {
2700 if(ftype<F_GB12 || ftype>F_GB14) {
2701
2702 if (interaction_function[ftype].flags & IF_BOND &&
2703 !(ftype == F_CONNBONDS || ftype == F_POSRES))
2704 {
2705 nbonds_np = idef->il[ftype].nr_nonperturbed;
2706 nbonds = idef->il[ftype].nr - nbonds_np;
2707 nat1 = interaction_function[ftype].nratoms + 1;
2708 if (nbonds > 0) {
2709 ind = interaction_function[ftype].nrnb_ind;
2710 iatom_fe = idef->il[ftype].iatoms + nbonds_np;
2711 dvdl = 0;
2712 if (ftype < F_LJ14 || ftype > F_LJC_PAIRS_NB) {
2713 v =
2714 interaction_function[ftype].ifunc(nbonds,iatom_fe,
2715 idef->iparams,
2716 (const rvec*)x,f,fr->fshift,
2717 pbc_null,g,lambda,&dvdl,md,fcd,
2718 global_atom_index);
2719 } else {
2720 v = do_listed_vdw_q(ftype,nbonds,iatom_fe,
2721 idef->iparams,
2722 (const rvec*)x,f,fr->fshift,
2723 pbc_null,g,
2724 lambda,&dvdl,
2725 md,fr,&enerd->grpp,global_atom_index);
2726 }
2727 if (ind != -1)
2728 inc_nrnb(nrnb,ind,nbonds/nat1);
2729 epot[ftype] += v;
2730 }
2731 }
2732 }
2733 }
2734
2735 sfree(fr->fshift);
2736 fr->fshift = fshift_orig;
2737 sfree(f);
2738 }