?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / gmxlib / bondfree.c
blobe8be08bc8d001c27326923add987a112c8b5716c
1 /*
2 * $Id$
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 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * Green Red Orange Magenta Azure Cyan Skyblue
28 */
29 static char *SRCID_bondfree_c = "$Id$";
30
31 #include <math.h>
32 #include "assert.h"
33 #include "physics.h"
34 #include "vec.h"
35 #include "vveclib.h"
36 #include "maths.h"
37 #include "txtdump.h"
38 #include "bondf.h"
39 #include "smalloc.h"
40 #include "pbc.h"
41 #include "ns.h"
42 #include "macros.h"
43 #include "names.h"
44 #include "fatal.h"
45 #include "mshift.h"
46 #include "main.h"
47 #include "disre.h"
48
49 static bool bPBC=FALSE;
50
51 void pbc_rvec_sub(rvec xi,rvec xj,rvec dx)
52 {
53 if (bPBC)
54 pbc_dx(xi,xj,dx);
55 else
56 rvec_sub(xi,xj,dx);
57 }
58
59 void calc_bonds(FILE *log,t_idef *idef,
60 rvec x_s[],rvec f[],
61 t_forcerec *fr,t_graph *g,
62 real epot[],t_nrnb *nrnb,
63 matrix box,real lambda,
64 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
65 {
66 static bool bFirst=TRUE;
67 int ftype,nbonds,ind,nat;
68
69 if (bFirst) {
70 bPBC = (getenv("GMXFULLPBC") != NULL);
71 if (bPBC)
72 fprintf(log,"Full PBC calculation = %s\n",bool_names[bPBC]);
73 bFirst = FALSE;
74 #ifdef DEBUG
75 p_graph(debug,"Bondage is fun",g);
76 #endif
77 }
78 for(ftype=0; (ftype<F_NRE); ftype++) {
79 if (interaction_function[ftype].flags & IF_BOND) {
80 nbonds=idef->il[ftype].nr;
81 if (nbonds > 0) {
82 epot[ftype]+=
83 interaction_function[ftype].ifunc(nbonds,idef->il[ftype].iatoms,
84 idef->iparams,x_s,f,fr,g,box,
85 lambda,&epot[F_DVDL],
86 md,ngrp,egnb,egcoul);
87 ind = interaction_function[ftype].nrnb_ind;
88 nat = interaction_function[ftype].nratoms+1;
89 if (ind != -1)
90 inc_nrnb(nrnb,ind,nbonds/nat);
91 }
92 }
93 }
94 }
95
96 /*
97 * Morse potential bond by Frank Everdij
98 *
99 * Three parameters needed:
100 *
101 * b0 = equilibrium distance in nm
102 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
103 * cb = well depth in kJ/mol
104 *
105 * Note: the potential is referenced to be +cb at infinite separation
106 * and zero at the equilibrium distance!
107 */
108
109 real morsebonds(int nbonds,
110 t_iatom forceatoms[],t_iparams forceparams[],
111 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
112 matrix box,real lambda,real *dvdl,
113 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
114 {
115 const real one=1.0;
116 const real two=2.0;
117 real dr,dr2,temp,omtemp,cbomtemp,fbond,vbond,fij,b0,be,cb,vtot;
118 rvec dx;
119 int i,m,ki,kj,type,ai,aj;
120
121 vtot = 0.0;
122 for(i=0; (i<nbonds); ) {
123 type = forceatoms[i++];
124 ai = forceatoms[i++];
125 aj = forceatoms[i++];
126
127 b0 = forceparams[type].morse.b0;
128 be = forceparams[type].morse.beta;
129 cb = forceparams[type].morse.cb;
130
131 pbc_rvec_sub(x[ai],x[aj],dx); /* 3 */
132 dr2 = iprod(dx,dx); /* 5 */
133 dr = sqrt(dr2); /* 10 */
134 temp = exp(-be*(dr-b0)); /* 12 */
135
136 if (temp == one)
137 continue;
138
139 omtemp = one-temp; /* 1 */
140 cbomtemp = cb*omtemp; /* 1 */
141 vbond = cbomtemp*omtemp; /* 1 */
142 fbond = -two*be*temp*cbomtemp*invsqrt(dr2); /* 9 */
143 vtot += vbond; /* 1 */
144
145 ki=SHIFT_INDEX(g,ai);
146 kj=SHIFT_INDEX(g,aj);
147 for (m=0; (m<DIM); m++) { /* 15 */
148 fij=fbond*dx[m];
149 f[ai][m]+=fij;
150 f[aj][m]-=fij;
151 fr->fshift[ki][m]+=fij;
152 fr->fshift[kj][m]-=fij;
153 }
154 } /* 58 TOTAL */
155 return vtot;
156 }
157
158 real harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
159 real *V,real *F)
160 {
161 const real half=0.5;
162 real L1,kk,x0,dx,dx2;
163 real v,f,dvdl;
164
165 L1 = 1.0-lambda;
166 kk = L1*kA+lambda*kB;
167 x0 = L1*xA+lambda*xB;
168
169 dx = x-x0;
170 dx2 = dx*dx;
171
172 f = -kk*dx;
173 v = half*kk*dx2;
174 dvdl = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
175
176 *F = f;
177 *V = v;
178
179 return dvdl;
180
181 /* That was 19 flops */
182 }
183
184
185 real bonds(int nbonds,
186 t_iatom forceatoms[],t_iparams forceparams[],
187 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
188 matrix box,real lambda,real *dvdlambda,
189 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
190 {
191 int i,m,ki,kj,ai,aj,type;
192 real dr,dr2,fbond,vbond,fij,vtot;
193 rvec dx;
194
195 vtot = 0.0;
196 for(i=0; (i<nbonds); ) {
197 type = forceatoms[i++];
198 ai = forceatoms[i++];
199 aj = forceatoms[i++];
200
201 pbc_rvec_sub(x[ai],x[aj],dx); /* 3 */
202 dr2=iprod(dx,dx); /* 5 */
203 dr=sqrt(dr2); /* 10 */
204
205 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
206 forceparams[type].harmonic.krB,
207 forceparams[type].harmonic.rA,
208 forceparams[type].harmonic.rB,
209 dr,lambda,&vbond,&fbond); /* 19 */
210
211 if (dr2 == 0.0)
212 continue;
213
214 vtot += vbond;/* 1*/
215 fbond *= invsqrt(dr2); /* 6 */
216 #ifdef DEBUG
217 if (debug)
218 fprintf(debug,"BONDS: dr = %10g vbond = %10g fbond = %10g\n",
219 dr,vbond,fbond);
220 #endif
221 ki=SHIFT_INDEX(g,ai);
222 kj=SHIFT_INDEX(g,aj);
223 for (m=0; (m<DIM); m++) { /* 15 */
224 fij=fbond*dx[m];
225 f[ai][m]+=fij;
226 f[aj][m]-=fij;
227 fr->fshift[ki][m]+=fij;
228 fr->fshift[kj][m]-=fij;
229 }
230 } /* 59 TOTAL */
231 return vtot;
232 }
233
234 real water_pol(int nbonds,
235 t_iatom forceatoms[],t_iparams forceparams[],
236 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
237 matrix box,real lambda,real *dvdlambda,
238 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
239 {
240 /* This routine implements anisotropic polarizibility for water, through
241 * a shell connected to a dummy with spring constant that differ in the
242 * three spatial dimensions in the molecular frame.
243 */
244 int i,m,aO,aH1,aH2,aD,aS,type;
245 rvec dOH1,dOH2,dHH,dOD,dDS,nW,kk,dx,kdx,proj;
246 #ifdef DEBUG
247 rvec df;
248 #endif
249 real vtot,fij,r_HH,r_OD,r_nW,tx,ty,tz;
250
251 vtot = 0.0;
252 if (nbonds > 0) {
253 type = forceatoms[0];
254 kk[XX] = forceparams[type].wpol.kx;
255 kk[YY] = forceparams[type].wpol.ky;
256 kk[ZZ] = forceparams[type].wpol.kz;
257 r_HH = 1.0/forceparams[type].wpol.rHH;
258 r_OD = 1.0/forceparams[type].wpol.rOD;
259 if (debug) {
260 fprintf(debug,"WPOL: kk = %10.3f %10.3f %10.3f\n",
261 kk[XX],kk[YY],kk[ZZ]);
262 fprintf(debug,"WPOL: rOH = %10.3f rHH = %10.3f rOD = %10.3f\n",
263 forceparams[type].wpol.rOH,
264 forceparams[type].wpol.rHH,
265 forceparams[type].wpol.rOD);
266 }
267 for(i=0; (i<nbonds); ) {
268 type = forceatoms[i++];
269 aO = forceatoms[i++];
270 aH1 = aO+1;
271 aH2 = aO+2;
272 aD = aO+3;
273 aS = aO+4;
274
275 /* Compute vectors describing the water frame */
276 rvec_sub(x[aH1],x[aO], dOH1);
277 rvec_sub(x[aH2],x[aO], dOH2);
278 rvec_sub(x[aH2],x[aH1],dHH);
279 rvec_sub(x[aD], x[aO], dOD);
280 rvec_sub(x[aS], x[aD], dDS);
281 oprod(dOH1,dOH2,nW);
282
283 /* Compute inverse length of normal vector
284 * (this one could be precomputed, but I'm too lazy now)
285 */
286 r_nW = invsqrt(iprod(nW,nW));
287 /* This is for precision, but does not make a big difference,
288 * it can go later.
289 */
290 r_OD = invsqrt(iprod(dOD,dOD));
291
292 /* Normalize the vectors in the water frame */
293 svmul(r_nW,nW,nW);
294 svmul(r_HH,dHH,dHH);
295 svmul(r_OD,dOD,dOD);
296
297 /* Compute displacement of shell along components of the vector */
298 dx[ZZ] = iprod(dDS,dOD);
299 /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
300 for(m=0; (m<DIM); m++)
301 proj[m] = dDS[m]-dx[ZZ]*dOD[m];
302
303 /*dx[XX] = iprod(dDS,nW);
304 dx[YY] = iprod(dDS,dHH);*/
305 dx[XX] = iprod(proj,nW);
306 for(m=0; (m<DIM); m++)
307 proj[m] -= dx[XX]*nW[m];
308 dx[YY] = iprod(proj,dHH);
309 /*#define DEBUG*/
310 #ifdef DEBUG
311 if (debug) {
312 fprintf(debug,"WPOL: dx2=%10g dy2=%10g dz2=%10g sum=%10g dDS^2=%10g\n",
313 sqr(dx[XX]),sqr(dx[YY]),sqr(dx[ZZ]),iprod(dx,dx),iprod(dDS,dDS));
314 fprintf(debug,"WPOL: dHH=(%10g,%10g,%10g)\n",dHH[XX],dHH[YY],dHH[ZZ]);
315 fprintf(debug,"WPOL: dOD=(%10g,%10g,%10g), 1/r_OD = %10g\n",
316 dOD[XX],dOD[YY],dOD[ZZ],1/r_OD);
317 fprintf(debug,"WPOL: nW =(%10g,%10g,%10g), 1/r_nW = %10g\n",
318 nW[XX],nW[YY],nW[ZZ],1/r_nW);
319 fprintf(debug,"WPOL: dx =%10g, dy =%10g, dz =%10g\n",
320 dx[XX],dx[YY],dx[ZZ]);
321 fprintf(debug,"WPOL: dDSx=%10g, dDSy=%10g, dDSz=%10g\n",
322 dDS[XX],dDS[YY],dDS[ZZ]);
323 }
324 #endif
325 /* Now compute the forces and energy */
326 kdx[XX] = kk[XX]*dx[XX];
327 kdx[YY] = kk[YY]*dx[YY];
328 kdx[ZZ] = kk[ZZ]*dx[ZZ];
329 vtot += iprod(dx,kdx);
330 for(m=0; (m<DIM); m++) {
331 /* This is a tensor operation but written out for speed */
332 tx = nW[m]*kdx[XX];
333 ty = dHH[m]*kdx[YY];
334 tz = dOD[m]*kdx[ZZ];
335 fij = -tx-ty-tz;
336 #ifdef DEBUG
337 df[m] = fij;
338 #endif
339 f[aS][m] += fij;
340 f[aD][m] -= fij;
341 }
342 #ifdef DEBUG
343 if (debug) {
344 fprintf(debug,"WPOL: vwpol=%g\n",0.5*iprod(dx,kdx));
345 fprintf(debug,"WPOL: df = (%10g, %10g, %10g)\n",df[XX],df[YY],df[ZZ]);
346 }
347 #endif
348 }
349 }
350 return 0.5*vtot;
351 }
352
353 real bond_angle(matrix box,
354 rvec xi,rvec xj,rvec xk, /* in */
355 rvec r_ij,rvec r_kj,real *costh) /* out */
356 /* Return value is the angle between the bonds i-j and j-k */
357 {
358 /* 41 FLOPS */
359 real th;
360
361 pbc_rvec_sub(xi,xj,r_ij); /* 3 */
362 pbc_rvec_sub(xk,xj,r_kj); /* 3 */
363
364 *costh=cos_angle(r_ij,r_kj); /* 25 */
365 th=acos(*costh); /* 10 */
366 /* 41 TOTAL */
367 return th;
368 }
369
370 real angles(int nbonds,
371 t_iatom forceatoms[],t_iparams forceparams[],
372 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
373 matrix box,real lambda,real *dvdlambda,
374 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
375 {
376 int i,ai,aj,ak,t,type;
377 rvec r_ij,r_kj;
378 real cos_theta,theta,dVdt,va,vtot;
379
380 vtot = 0.0;
381 for(i=0; (i<nbonds); ) {
382 type = forceatoms[i++];
383 ai = forceatoms[i++];
384 aj = forceatoms[i++];
385 ak = forceatoms[i++];
386
387 theta = bond_angle(box,x[ai],x[aj],x[ak],
388 r_ij,r_kj,&cos_theta); /* 41 */
389
390 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
391 forceparams[type].harmonic.krB,
392 forceparams[type].harmonic.rA*DEG2RAD,
393 forceparams[type].harmonic.rB*DEG2RAD,
394 theta,lambda,&va,&dVdt); /* 21 */
395 vtot += va;
396
397 {
398 int m;
399 real snt,st,sth;
400 real cik,cii,ckk;
401 real nrkj2,nrij2;
402 rvec f_i,f_j,f_k;
403
404 snt=sin(theta); /* 10 */
405 if (fabs(snt) < 1e-12)
406 snt=1e-12;
407 st = dVdt/snt; /* 10 */
408 sth = st*cos_theta; /* 1 */
409 #ifdef DEBUG
410 if (debug)
411 fprintf(debug,"ANGLES: theta = %10g vth = %10g dV/dtheta = %10g\n",
412 theta*RAD2DEG,va,dVdt);
413 #endif
414 nrkj2=iprod(r_kj,r_kj); /* 5 */
415 nrij2=iprod(r_ij,r_ij);
416
417 cik=st*invsqrt(nrkj2*nrij2); /* 12 */
418 cii=sth/nrij2; /* 10 */
419 ckk=sth/nrkj2; /* 10 */
420
421 for (m=0; (m<DIM); m++) { /* 39 */
422 f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
423 f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
424 f_j[m]=-f_i[m]-f_k[m];
425 f[ai][m]+=f_i[m];
426 f[aj][m]+=f_j[m];
427 f[ak][m]+=f_k[m];
428 }
429 t=SHIFT_INDEX(g,ai);
430 rvec_inc(fr->fshift[t],f_i);
431 t=SHIFT_INDEX(g,aj);
432 rvec_inc(fr->fshift[t],f_j);
433 t=SHIFT_INDEX(g,ak);
434 rvec_inc(fr->fshift[t],f_k);
435 } /* 168 TOTAL */
436 }
437 return vtot;
438 }
439
440 real dih_angle(matrix box,
441 rvec xi,rvec xj,rvec xk,rvec xl,
442 rvec r_ij,rvec r_kj,rvec r_kl,rvec m,rvec n,
443 real *cos_phi,real *sign)
444 {
445 real ipr,phi;
446
447 pbc_rvec_sub(xi,xj,r_ij); /* 3 */
448 pbc_rvec_sub(xk,xj,r_kj); /* 3 */
449 pbc_rvec_sub(xk,xl,r_kl); /* 3 */
450
451 oprod(r_ij,r_kj,m); /* 9 */
452 oprod(r_kj,r_kl,n); /* 9 */
453 *cos_phi=cos_angle(m,n); /* 41 */
454 phi=acos(*cos_phi); /* 10 */
455 ipr=iprod(r_ij,n); /* 5 */
456 (*sign)=(ipr<0.0)?-1.0:1.0;
457 phi=(*sign)*phi; /* 1 */
458 /* 84 TOTAL */
459 return phi;
460 }
461
462 void do_dih_fup(int i,int j,int k,int l,real ddphi,
463 rvec r_ij,rvec r_kj,rvec r_kl,
464 rvec m,rvec n,rvec f[],t_forcerec *fr,t_graph *g,
465 rvec x[])
466 {
467 /* 143 FLOPS */
468 rvec f_i,f_j,f_k,f_l;
469 rvec uvec,vvec,svec;
470 real ipr,nrkj,nrkj2;
471 real a,p,q;
472 int t;
473
474 ipr = iprod(m,m); /* 5 */
475 nrkj2 = iprod(r_kj,r_kj); /* 5 */
476 nrkj = sqrt(nrkj2); /* 10 */
477 a = -ddphi*nrkj/ipr; /* 11 */
478 svmul(a,m,f_i); /* 3 */
479 ipr = iprod(n,n); /* 5 */
480 a = ddphi*nrkj/ipr; /* 11 */
481 svmul(a,n,f_l); /* 3 */
482 p = iprod(r_ij,r_kj); /* 5 */
483 p /= nrkj2; /* 10 */
484 q = iprod(r_kl,r_kj); /* 5 */
485 q /= nrkj2; /* 10 */
486 svmul(p,f_i,uvec); /* 3 */
487 svmul(q,f_l,vvec); /* 3 */
488 rvec_sub(uvec,vvec,svec); /* 3 */
489 rvec_sub(f_i,svec,f_j); /* 3 */
490 rvec_add(f_l,svec,f_k); /* 3 */
491 rvec_inc(f[i],f_i); /* 3 */
492 rvec_dec(f[j],f_j); /* 3 */
493 rvec_dec(f[k],f_k); /* 3 */
494 rvec_inc(f[l],f_l); /* 3 */
495
496 t=SHIFT_INDEX(g,i);
497 rvec_inc(fr->fshift[t],f_i);
498 t=SHIFT_INDEX(g,j);
499 rvec_dec(fr->fshift[t],f_j);
500 t=SHIFT_INDEX(g,k);
501 rvec_dec(fr->fshift[t],f_k);
502 t=SHIFT_INDEX(g,l);
503 rvec_inc(fr->fshift[t],f_l);
504 /* 112 TOTAL */
505 }
506
507
508 real dopdihs(real cpA,real cpB,real phiA,real phiB,int mult,
509 real phi,real lambda,real *V,real *F)
510 {
511 real v,dvdl,mdphi,v1,sdphi,ddphi;
512 real L1 = 1.0-lambda;
513 real ph0 = DEG2RAD*(L1*phiA+lambda*phiB);
514 real cp = L1*cpA + lambda*cpB;
515
516 mdphi = mult*phi-ph0;
517 sdphi = sin(mdphi);
518 ddphi = -cp*mult*sdphi;
519 v1 = 1.0+cos(mdphi);
520 v = cp*v1;
521
522 dvdl = (cpB-cpA)*v1 - cp*(phiA-phiB)*sdphi;
523
524 *V = v;
525 *F = ddphi;
526
527 return dvdl;
528
529 /* That was 40 flops */
530 }
531
532 static real dopdihs_min(real cpA,real cpB,real phiA,real phiB,int mult,
533 real phi,real lambda,real *V,real *F)
534 /* similar to dopdihs, except for a minus sign *
535 * and a different treatment of mult/phi0 */
536 {
537 real v,dvdl,mdphi,v1,sdphi,ddphi;
538 real L1 = 1.0-lambda;
539 real ph0 = DEG2RAD*(L1*phiA+lambda*phiB);
540 real cp = L1*cpA + lambda*cpB;
541
542 mdphi = mult*(phi-ph0);
543 sdphi = sin(mdphi);
544 ddphi = cp*mult*sdphi;
545 v1 = 1.0-cos(mdphi);
546 v = cp*v1;
547
548 dvdl = (cpB-cpA)*v1 - cp*(phiA-phiB)*sdphi;
549
550 *V = v;
551 *F = ddphi;
552
553 return dvdl;
554
555 /* That was 40 flops */
556 }
557
558 real pdihs(int nbonds,
559 t_iatom forceatoms[],t_iparams forceparams[],
560 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
561 matrix box,real lambda,real *dvdlambda,
562 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
563 {
564 int i,type,ai,aj,ak,al;
565 rvec r_ij,r_kj,r_kl,m,n;
566 real phi,cos_phi,sign,ddphi,vpd,vtot;
567
568 vtot = 0.0;
569 for(i=0; (i<nbonds); ) {
570 type = forceatoms[i++];
571 ai = forceatoms[i++];
572 aj = forceatoms[i++];
573 ak = forceatoms[i++];
574 al = forceatoms[i++];
575
576 phi=dih_angle(box,x[ai],x[aj],x[ak],x[al],r_ij,r_kj,r_kl,m,n,
577 &cos_phi,&sign); /* 84 */
578
579 *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
580 forceparams[type].pdihs.cpB,
581 forceparams[type].pdihs.phiA,
582 forceparams[type].pdihs.phiB,
583 forceparams[type].pdihs.mult,
584 phi,lambda,&vpd,&ddphi);
585
586 vtot += vpd;
587 do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
588 f,fr,g,x); /* 112 */
589
590 #ifdef DEBUG
591 fprintf(debug,"pdih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
592 ai,aj,ak,al,cos_phi,phi);
593 #endif
594 } /* 223 TOTAL */
595
596 return vtot;
597 }
598
599 real idihs(int nbonds,
600 t_iatom forceatoms[],t_iparams forceparams[],
601 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
602 matrix box,real lambda,real *dvdlambda,
603 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
604 {
605 int i,type,ai,aj,ak,al;
606 real phi,cos_phi,ddphi,sign,vid,vtot;
607 rvec r_ij,r_kj,r_kl,m,n;
608
609 vtot = 0.0;
610 for(i=0; (i<nbonds); ) {
611 type = forceatoms[i++];
612 ai = forceatoms[i++];
613 aj = forceatoms[i++];
614 ak = forceatoms[i++];
615 al = forceatoms[i++];
616
617 phi=dih_angle(box,x[ai],x[aj],x[ak],x[al],r_ij,r_kj,r_kl,m,n,
618 &cos_phi,&sign); /* 84 */
619
620 *dvdlambda += harmonic(forceparams[type].harmonic.krA,
621 forceparams[type].harmonic.krB,
622 forceparams[type].harmonic.rA*DEG2RAD,
623 forceparams[type].harmonic.rB*DEG2RAD,
624 phi,lambda,&vid,&ddphi); /* 21 */
625
626 vtot += vid;
627 do_dih_fup(ai,aj,ak,al,(real)(-ddphi),r_ij,r_kj,r_kl,m,n,
628 f,fr,g,x); /* 112 */
629 /* 217 TOTAL */
630 #ifdef DEBUG
631 fprintf("idih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
632 ai,aj,ak,al,cos_phi,phi);
633 #endif
634 }
635 return vtot;
636 }
637
638 real posres(int nbonds,
639 t_iatom forceatoms[],t_iparams forceparams[],
640 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
641 matrix box,real lambda,real *dvdlambda,
642 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
643 {
644 int i,ai,m,type;
645 real v,vtot,fi,*fc;
646 rvec dx;
647
648 vtot = 0.0;
649 for(i=0; (i<nbonds); ) {
650 type = forceatoms[i++];
651 ai = forceatoms[i++];
652 fc = forceparams[type].posres.fc;
653
654 pbc_dx(x[ai],forceparams[type].posres.pos0,dx);
655 v=0;
656 for (m=0; (m<DIM); m++) {
657 fi = f[ai][m] - fc[m]*dx[m];
658 v += 0.5*fc[m]*dx[m]*dx[m];
659 f[ai][m] = fi;
660 }
661 vtot += v;
662 }
663 return vtot;
664 }
665
666 static real low_angres(int nbonds,
667 t_iatom forceatoms[],t_iparams forceparams[],
668 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
669 matrix box,real lambda,real *dvdlambda,
670 bool bZAxis)
671 {
672 int i,m,type,ai,aj,ak,al,t;
673 real phi,cos_phi,sign,vid,vtot,dVdphi;
674 rvec r_ij,r_kl,f_i,f_k;
675 real st,sth,sin_phi,nrij2,nrkl2,c,cij,ckl;
676
677 vtot = 0.0;
678 ak=al=0; /* to avoid warnings */
679 for(i=0; i<nbonds; ) {
680 type = forceatoms[i++];
681 ai = forceatoms[i++];
682 aj = forceatoms[i++];
683 pbc_rvec_sub(x[aj],x[ai],r_ij); /* 3 */
684 if (!bZAxis) {
685 ak = forceatoms[i++];
686 al = forceatoms[i++];
687 pbc_rvec_sub(x[al],x[ak],r_kl); /* 3 */
688 } else {
689 r_kl[XX] = 0;
690 r_kl[YY] = 0;
691 r_kl[ZZ] = 1;
692 }
693
694 cos_phi = cos_angle(r_ij,r_kl); /* 25 */
695 phi = acos(cos_phi); /* 10 */
696
697 *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
698 forceparams[type].pdihs.cpB,
699 forceparams[type].pdihs.phiA,
700 forceparams[type].pdihs.phiB,
701 forceparams[type].pdihs.mult,
702 phi,lambda,&vid,&dVdphi); /* 40 */
703
704 vtot += vid;
705
706 sin_phi = sin(phi); /* 10 */
707 if (fabs(sin_phi) < 1e-12)
708 sin_phi=1e-12;
709 st = -dVdphi/sin_phi; /* 10 */
710 sth = st*cos_phi; /* 1 */
711 nrij2 = iprod(r_ij,r_ij); /* 5 */
712 nrkl2 = iprod(r_kl,r_kl); /* 5 */
713
714 c = st*invsqrt(nrij2*nrkl2); /* 11 */
715 cij = sth/nrij2; /* 10 */
716 ckl = sth/nrkl2; /* 10 */
717
718 for (m=0; m<DIM; m++) { /* 18+18 */
719 f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
720 f[ai][m] += f_i[m];
721 f[aj][m] -= f_i[m];
722 if (!bZAxis) {
723 f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
724 f[ak][m] += f_k[m];
725 f[al][m] -= f_k[m];
726 }
727 }
728 t=SHIFT_INDEX(g,ai);
729 rvec_inc(fr->fshift[t],f_i);
730 t=SHIFT_INDEX(g,aj);
731 rvec_dec(fr->fshift[t],f_i);
732 if (!bZAxis) {
733 t=SHIFT_INDEX(g,ak);
734 rvec_inc(fr->fshift[t],f_k);
735 t=SHIFT_INDEX(g,al);
736 rvec_dec(fr->fshift[t],f_k);
737 }
738 }
739
740 return vtot; /* 191 / 164 (bZAxis) total */
741 }
742
743 real angres(int nbonds,
744 t_iatom forceatoms[],t_iparams forceparams[],
745 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
746 matrix box,real lambda,real *dvdlambda,
747 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
748 {
749 return low_angres(nbonds,forceatoms,forceparams,x,f,fr,g,box,
750 lambda,dvdlambda,FALSE);
751 }
752
753 real angresz(int nbonds,
754 t_iatom forceatoms[],t_iparams forceparams[],
755 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
756 matrix box,real lambda,real *dvdlambda,
757 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
758 {
759 return low_angres(nbonds,forceatoms,forceparams,x,f,fr,g,box,
760 lambda,dvdlambda,TRUE);
761 }
762
763 real unimplemented(int nbonds,
764 t_iatom forceatoms[],t_iparams forceparams[],
765 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
766 matrix box,real lambda,real *dvdlambda,
767 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
768 {
769 fatal_error(0,"*** you are using a not implemented function");
770
771 return 0.0; /* To make the compiler happy */
772 }
773
774 static void my_fatal(char *s,int ai,int aj,int ak,int al)
775 {
776 fatal_error(0,"%s is NaN in rbdih, ai-ak=%d,%d,%d,%d",s,ai,aj,ak,al);
777 }
778
779 #define CHECK(s) if ((s != s) || (s < -1e10) || (s > 1e10)) my_fatal(#s,ai,aj,ak,al)
780
781 real rbdihs(int nbonds,
782 t_iatom forceatoms[],t_iparams forceparams[],
783 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
784 matrix box,real lambda,real *dvdlambda,
785 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
786 {
787 static const real c0=0.0,c1=1.0,c2=2.0,c3=3.0,c4=4.0,c5=5.0;
788 int type,ai,aj,ak,al,i,j;
789 rvec r_ij,r_kj,r_kl,m,n;
790 real parm[NR_RBDIHS];
791 real phi,cos_phi,rbp;
792 real v,sign,ddphi,sin_phi;
793 real cosfac,vtot;
794
795 vtot = 0.0;
796 for(i=0; (i<nbonds); ) {
797 type = forceatoms[i++];
798 ai = forceatoms[i++];
799 aj = forceatoms[i++];
800 ak = forceatoms[i++];
801 al = forceatoms[i++];
802
803 phi=dih_angle(box,x[ai],x[aj],x[ak],x[al],r_ij,r_kj,r_kl,m,n,
804 &cos_phi,&sign); /* 84 */
805
806 /* Change to polymer convention */
807 if (phi < c0)
808 phi += M_PI;
809 else
810 phi -= M_PI; /* 1 */
811 cos_phi = -cos_phi; /* 1 */
812
813 sin_phi=sin(phi);
814
815 for(j=0; (j<NR_RBDIHS); j++)
816 parm[j] = forceparams[type].rbdihs.rbc[j];
817
818 /* Calculate cosine powers */
819 /* Calculate the energy */
820 /* Calculate the derivative */
821 v = parm[0];
822 ddphi = c0;
823 cosfac = c1;
824
825 rbp = parm[1];
826 ddphi += rbp*cosfac;
827 cosfac *= cos_phi;
828 v += cosfac*rbp;
829 rbp = parm[2];
830 ddphi += c2*rbp*cosfac;
831 cosfac *= cos_phi;
832 v += cosfac*rbp;
833 rbp = parm[3];
834 ddphi += c3*rbp*cosfac;
835 cosfac *= cos_phi;
836 v += cosfac*rbp;
837 rbp = parm[4];
838 ddphi += c4*rbp*cosfac;
839 cosfac *= cos_phi;
840 v += cosfac*rbp;
841 rbp = parm[5];
842 ddphi += c5*rbp*cosfac;
843 cosfac *= cos_phi;
844 v += cosfac*rbp;
845
846 ddphi = -ddphi*sin_phi; /* 11 */
847
848 do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
849 f,fr,g,x); /* 112 */
850 vtot += v;
851 }
852 return vtot;
853 }
854
855 /***********************************************************
856 *
857 * G R O M O S 9 6 F U N C T I O N S
858 *
859 ***********************************************************/
860 real g96harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
861 real *V,real *F)
862 {
863 const real half=0.5;
864 real L1,kk,x0,dx,dx2;
865 real v,f,dvdl;
866
867 L1 = 1.0-lambda;
868 kk = L1*kA+lambda*kB;
869 x0 = L1*xA+lambda*xB;
870
871 dx = x-x0;
872 dx2 = dx*dx;
873
874 f = -kk*dx;
875 v = half*kk*dx2;
876 dvdl = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
877
878 *F = f;
879 *V = v;
880
881 return dvdl;
882
883 /* That was 21 flops */
884 }
885
886 real g96bonds(int nbonds,
887 t_iatom forceatoms[],t_iparams forceparams[],
888 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
889 matrix box,real lambda,real *dvdlambda,
890 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
891 {
892 int i,m,ki,kj,ai,aj,type;
893 real dr2,fbond,vbond,fij,vtot;
894 rvec dx;
895
896 vtot = 0.0;
897 for(i=0; (i<nbonds); ) {
898 type = forceatoms[i++];
899 ai = forceatoms[i++];
900 aj = forceatoms[i++];
901
902 pbc_rvec_sub(x[ai],x[aj],dx); /* 3 */
903 dr2=iprod(dx,dx); /* 5 */
904
905 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
906 forceparams[type].harmonic.krB,
907 forceparams[type].harmonic.rA,
908 forceparams[type].harmonic.rB,
909 dr2,lambda,&vbond,&fbond);
910
911 vtot += 0.5*vbond; /* 1*/
912 #ifdef DEBUG
913 if (debug)
914 fprintf(debug,"G96-BONDS: dr = %10g vbond = %10g fbond = %10g\n",
915 sqrt(dr2),vbond,fbond);
916 #endif
917 ki=SHIFT_INDEX(g,ai);
918 kj=SHIFT_INDEX(g,aj);
919 for (m=0; (m<DIM); m++) { /* 15 */
920 fij=fbond*dx[m];
921 f[ai][m]+=fij;
922 f[aj][m]-=fij;
923 fr->fshift[ki][m]+=fij;
924 fr->fshift[kj][m]-=fij;
925 }
926 } /* 44 TOTAL */
927 return vtot;
928 }
929
930 real g96bond_angle(matrix box,
931 rvec xi,rvec xj,rvec xk, /* in */
932 rvec r_ij,rvec r_kj) /* out */
933 /* Return value is the angle between the bonds i-j and j-k */
934 {
935 real costh;
936
937 pbc_rvec_sub(xi,xj,r_ij); /* 3 */
938 pbc_rvec_sub(xk,xj,r_kj); /* 3 */
939
940 costh=cos_angle(r_ij,r_kj); /* 25 */
941 /* 41 TOTAL */
942 return costh;
943 }
944
945 real g96angles(int nbonds,
946 t_iatom forceatoms[],t_iparams forceparams[],
947 rvec x[],rvec f[],t_forcerec *fr,t_graph *g,
948 matrix box,real lambda,real *dvdlambda,
949 t_mdatoms *md,int ngrp,real egnb[],real egcoul[])
950 {
951 int i,ai,aj,ak,t,type,m;
952 rvec r_ij,r_kj;
953 real cos_theta,dVdt,va,vtot;
954 real rij_1,rij_2,rkj_1,rkj_2,rijrkj_1;
955 rvec f_i,f_j,f_k;
956
957 vtot = 0.0;
958 for(i=0; (i<nbonds); ) {
959 type = forceatoms[i++];
960 ai = forceatoms[i++];
961 aj = forceatoms[i++];
962 ak = forceatoms[i++];
963
964 cos_theta = g96bond_angle(box,x[ai],x[aj],x[ak],r_ij,r_kj);
965
966 *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
967 forceparams[type].harmonic.krB,
968 forceparams[type].harmonic.rA,
969 forceparams[type].harmonic.rB,
970 cos_theta,lambda,&va,&dVdt);
971 vtot += va;
972
973 rij_1 = invsqrt(iprod(r_ij,r_ij));
974 rkj_1 = invsqrt(iprod(r_kj,r_kj));
975 rij_2 = rij_1*rij_1;
976 rkj_2 = rkj_1*rkj_1;
977 rijrkj_1 = rij_1*rkj_1; /* 23 */
978
979 #ifdef DEBUG
980 if (debug)
981 fprintf(debug,"G96ANGLES: costheta = %10g vth = %10g dV/dct = %10g\n",
982 cos_theta,va,dVdt);
983 #endif
984 for (m=0; (m<DIM); m++) { /* 42 */
985 f_i[m]=dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
986 f_k[m]=dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
987 f_j[m]=-f_i[m]-f_k[m];
988 f[ai][m]+=f_i[m];
989 f[aj][m]+=f_j[m];
990 f[ak][m]+=f_k[m];
991 }
992 t = SHIFT_INDEX(g,ai);
993 rvec_inc(fr->fshift[t],f_i);
994 t = SHIFT_INDEX(g,aj);
995 rvec_inc(fr->fshift[t],f_j);
996 t = SHIFT_INDEX(g,ak);
997 rvec_inc(fr->fshift[t],f_k); /* 9 */
998 /* 163 TOTAL */
999 }
1000 return vtot;
1001 }
1002
got merge a long time ago