+real linear_angles(int nbonds,
+ const t_iatom forceatoms[],const t_iparams forceparams[],
+ const rvec x[],rvec f[],rvec fshift[],
+ const t_pbc *pbc,const t_graph *g,
+ real lambda,real *dvdlambda,
+ const t_mdatoms *md,t_fcdata *fcd,
+ int *global_atom_index)
+{
+ int i,m,ai,aj,ak,t1,t2,type;
+ rvec f_i,f_j,f_k;
+ real L1,kA,kB,aA,aB,dr,dr2,va,vtot,a,b,klin,dvdl;
+ ivec jt,dt_ij,dt_kj;
+ rvec r_ij,r_kj,r_ik,dx;
+
+ L1 = 1-lambda;
+ vtot = 0.0;
+ dvdl = 0.0;
+ for(i=0; (i<nbonds); ) {
+ type = forceatoms[i++];
+ ai = forceatoms[i++];
+ aj = forceatoms[i++];
+ ak = forceatoms[i++];
+
+ kA = forceparams[type].linangle.klinA;
+ kB = forceparams[type].linangle.klinB;
+ klin = L1*kA + lambda*kB;
+
+ aA = forceparams[type].linangle.aA;
+ aB = forceparams[type].linangle.aB;
+ a = L1*aA+lambda*aB;
+ b = 1-a;
+
+ t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
+ t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
+ rvec_sub(r_ij,r_kj,r_ik);
+
+ dr2 = 0;
+ for(m=0; (m<DIM); m++)
+ {
+ dr = - a * r_ij[m] - b * r_kj[m];
+ dr2 += dr*dr;
+ dx[m] = dr;
+ f_i[m] = a*klin*dr;
+ f_k[m] = b*klin*dr;
+ f_j[m] = -(f_i[m]+f_k[m]);
+ f[ai][m] += f_i[m];
+ f[aj][m] += f_j[m];
+ f[ak][m] += f_k[m];
+ }
+ va = 0.5*klin*dr2;
+ dvdl += 0.5*(kB-kA)*dr2 + klin*(aB-aA)*iprod(dx,r_ik);
+
+ vtot += va;
+
+ if (g) {
+ copy_ivec(SHIFT_IVEC(g,aj),jt);
+
+ ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
+ ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
+ t1=IVEC2IS(dt_ij);
+ t2=IVEC2IS(dt_kj);
+ }
+ rvec_inc(fshift[t1],f_i);
+ rvec_inc(fshift[CENTRAL],f_j);
+ rvec_inc(fshift[t2],f_k);
+ } /* 57 TOTAL */
+ *dvdlambda = dvdl;
+ return vtot;
+}
+