Project

General

Profile

bondfree.c

gmx/src/gmxlib/bondfree.c - David van der Spoel, 11/19/2007 07:19 PM

 
1
/*
2
 * $Id: bondfree.c,v 1.85.2.10 2007/11/19 18:16:58 spoel Exp $
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.3.2
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-2007, 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 Machine for Chemical Simulation
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
static int pbc_rvec_sub(const t_pbc *pbc,const rvec xi,const rvec xj,rvec dx)
61
{
62
  if (pbc) {
63
    return pbc_dx(pbc,xi,xj,dx);
64
  }
65
  else {
66
    rvec_sub(xi,xj,dx);
67
    return CENTRAL;
68
  }
69
}
70

    
71
void calc_bonds(FILE *fplog,const t_commrec *cr,const t_commrec *mcr,
72
                const t_idef *idef,
73
                rvec x[],rvec f[],
74
                t_forcerec *fr,const t_pbc *pbc,const t_graph *g,
75
                real epot[],t_nrnb *nrnb,
76
                real lambda,
77
                const t_mdatoms *md,int ngrp,real egnb[],real egcoul[],
78
                t_fcdata *fcd,
79
                int step,bool bSepDVDL)
80
{
81
  int    ftype,nbonds,ind,nat;
82
  real   v,dvdl;
83
  const  t_pbc *pbc_null;
84

    
85
  if (fr->ePBC == epbcFULL)
86
    pbc_null = pbc;
87
  else
88
    pbc_null = NULL;
89

    
90
  if (bSepDVDL)
91
    fprintf(fplog,"Step %d: bonded V and dVdl for node %d:\n",step,cr->nodeid);
92

    
93
#ifdef DEBUG
94
  if (g && debug)
95
    p_graph(debug,"Bondage is fun",g);
96
#endif
97
  
98
  /* Do pre force calculation stuff which might require communication */
99
  if (idef->il[F_ORIRES].nr)
100
    epot[F_ORIRESDEV] = calc_orires_dev(mcr,idef->il[F_ORIRES].nr,
101
                                        idef->il[F_ORIRES].iatoms,
102
                                        idef->iparams,md,(const rvec*)x,
103
                                        pbc_null,fcd);
104
  if (idef->il[F_DISRES].nr)
105
    calc_disres_R_6(mcr,idef->il[F_DISRES].nr,
106
                    idef->il[F_DISRES].iatoms,
107
                    idef->iparams,(const rvec*)x,pbc_null,fcd);
108
  
109
  /* Loop over all bonded force types to calculate the bonded forces */
110
  for(ftype=0; (ftype<F_NRE); ftype++) {
111
    if (interaction_function[ftype].flags & IF_BOND && ftype!=F_CONNBONDS) {
112
      nbonds=idef->il[ftype].nr;
113
      if (nbonds > 0) {
114
        ind = interaction_function[ftype].nrnb_ind;
115
        nat = interaction_function[ftype].nratoms+1;
116
        dvdl = 0;
117
        if (ftype != F_LJ14) {
118
          v = interaction_function[ftype].ifunc
119
            (nbonds,idef->il[ftype].iatoms,
120
             idef->iparams,
121
             (const rvec*)x,f,fr->fshift,
122
             (ftype == F_POSRES) ? pbc : pbc_null,
123
             g,lambda,&dvdl,md,fcd);
124
          if (bSepDVDL) {
125
            fprintf(fplog,"  %-23s #%4d  V %12.5e  dVdl %12.5e\n",
126
                    interaction_function[ftype].longname,nbonds/nat,v,dvdl);
127
          }
128
        } else {
129
          v = do_nonbonded14(nbonds,idef->il[ftype].iatoms,
130
                             idef->iparams,
131
                             (const rvec*)x,f,fr->fshift,
132
                             pbc_null,g,
133
                             lambda,&dvdl,
134
                             md,fr,ngrp,egnb,egcoul);
135
          if (bSepDVDL) {
136
            fprintf(fplog,"  %-5s + %-15s #%4d                  dVdl %12.5e\n",
137
                    interaction_function[ftype].longname,
138
                    interaction_function[F_COUL14].longname,nbonds/nat,dvdl);
139
          }
140
        }
141
        if (ind != -1)
142
          inc_nrnb(nrnb,ind,nbonds/nat);
143
        epot[ftype]  += v;
144
        epot[F_DVDL] += dvdl;
145
      }
146
    }
147
  }
148
  /* Copy the sum of violations for the distance restraints from fcd */
149
  epot[F_DISRESVIOL] = fcd->disres.sumviol;
150
}
151

    
152
/*
153
 * Morse potential bond by Frank Everdij
154
 *
155
 * Three parameters needed:
156
 *
157
 * b0 = equilibrium distance in nm
158
 * be = beta in nm^-1 (actually, it's nu_e*Sqrt(2*pi*pi*mu/D_e))
159
 * cb = well depth in kJ/mol
160
 *
161
 * Note: the potential is referenced to be +cb at infinite separation
162
 *       and zero at the equilibrium distance!
163
 */
164

    
165
real morse_bonds(int nbonds,
166
                 const t_iatom forceatoms[],const t_iparams forceparams[],
167
                 const rvec x[],rvec f[],rvec fshift[],
168
                 const t_pbc *pbc,const t_graph *g,
169
                 real lambda,real *dvdl,
170
                 const t_mdatoms *md,t_fcdata *fcd)
171
{
172
  const real one=1.0;
173
  const real two=2.0;
174
  real  dr,dr2,temp,omtemp,cbomtemp,fbond,vbond,fij,b0,be,cb,vtot;
175
  rvec  dx;
176
  int   i,m,ki,type,ai,aj;
177
  ivec  dt;
178

    
179
  vtot = 0.0;
180
  for(i=0; (i<nbonds); ) {
181
    type = forceatoms[i++];
182
    ai   = forceatoms[i++];
183
    aj   = forceatoms[i++];
184
    
185
    b0   = forceparams[type].morse.b0;
186
    be   = forceparams[type].morse.beta;
187
    cb   = forceparams[type].morse.cb;
188

    
189
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
190
    dr2  = iprod(dx,dx);                            /*   5          */
191
    dr   = dr2*invsqrt(dr2);                        /*  10          */
192
    temp = exp(-be*(dr-b0));                        /*  12          */
193
    
194
    if (temp == one)
195
      continue;
196

    
197
    omtemp   = one-temp;                               /*   1          */
198
    cbomtemp = cb*omtemp;                              /*   1          */
199
    vbond    = cbomtemp*omtemp;                        /*   1          */
200
    fbond    = -two*be*temp*cbomtemp*invsqrt(dr2);      /*   9          */
201
    vtot    += vbond;       /* 1 */
202
    
203
    if (g) {
204
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
205
      ki = IVEC2IS(dt);
206
    }
207

    
208
    for (m=0; (m<DIM); m++) {                          /*  15          */
209
      fij=fbond*dx[m];
210
      f[ai][m]+=fij;
211
      f[aj][m]-=fij;
212
      fshift[ki][m]+=fij;
213
      fshift[CENTRAL][m]-=fij;
214
    }
215
  }                                           /*  58 TOTAL    */
216
  return vtot;
217
}
218

    
219
real cubic_bonds(int nbonds,
220
                 const t_iatom forceatoms[],const t_iparams forceparams[],
221
                 const rvec x[],rvec f[],rvec fshift[],
222
                 const t_pbc *pbc,const t_graph *g,
223
                 real lambda,real *dvdl,
224
                 const t_mdatoms *md,t_fcdata *fcd)
225
{
226
  const real three = 3.0;
227
  const real two   = 2.0;
228
  real  kb,b0,kcub;
229
  real  dr,dr2,dist,kdist,kdist2,fbond,vbond,fij,vtot;
230
  rvec  dx;
231
  int   i,m,ki,type,ai,aj;
232
  ivec  dt;
233

    
234
  vtot = 0.0;
235
  for(i=0; (i<nbonds); ) {
236
    type = forceatoms[i++];
237
    ai   = forceatoms[i++];
238
    aj   = forceatoms[i++];
239
    
240
    b0   = forceparams[type].cubic.b0;
241
    kb   = forceparams[type].cubic.kb;
242
    kcub = forceparams[type].cubic.kcub;
243

    
244
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);                /*   3          */
245
    dr2  = iprod(dx,dx);                                /*   5          */
246
    
247
    if (dr2 == 0.0)
248
      continue;
249
      
250
    dr         = dr2*invsqrt(dr2);                      /*  10          */
251
    dist       = dr-b0;
252
    kdist      = kb*dist;
253
    kdist2     = kdist*dist;
254
    
255
    vbond      = kdist2 + kcub*kdist2*dist;
256
    fbond      = -(two*kdist + three*kdist2*kcub)/dr;
257

    
258
    vtot      += vbond;       /* 21 */
259
    
260
    if (g) {
261
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
262
      ki=IVEC2IS(dt);
263
    }
264
    for (m=0; (m<DIM); m++) {                          /*  15          */
265
      fij=fbond*dx[m];
266
      f[ai][m]+=fij;
267
      f[aj][m]-=fij;
268
      fshift[ki][m]+=fij;
269
      fshift[CENTRAL][m]-=fij;
270
    }
271
  }                                           /*  54 TOTAL    */
272
  return vtot;
273
}
274

    
275
real FENE_bonds(int nbonds,
276
                const t_iatom forceatoms[],const t_iparams forceparams[],
277
                const rvec x[],rvec f[],rvec fshift[],
278
                const t_pbc *pbc,const t_graph *g,
279
                real lambda,real *dvdl,
280
                const t_mdatoms *md,t_fcdata *fcd)
281
{
282
  const real half=0.5;
283
  const real one=1.0;
284
  real  bm,kb;
285
  real  dr,dr2,bm2,omdr2obm2,fbond,vbond,fij,vtot;
286
  rvec  dx;
287
  int   i,m,ki,type,ai,aj;
288
  ivec  dt;
289

    
290
  vtot = 0.0;
291
  for(i=0; (i<nbonds); ) {
292
    type = forceatoms[i++];
293
    ai   = forceatoms[i++];
294
    aj   = forceatoms[i++];
295
    
296
    bm   = forceparams[type].fene.bm;
297
    kb   = forceparams[type].fene.kb;
298

    
299
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);            /*   3          */
300
    dr2  = iprod(dx,dx);                                /*   5          */
301
    
302
    if (dr2 == 0.0)
303
      continue;
304

    
305
    bm2 = bm*bm;
306

    
307
    if (dr2 >= bm2)
308
      gmx_fatal(FARGS,
309
                "r^2 (%f) >= bm^2 (%f) in FENE bond between atoms %d and %d",
310
                dr2,bm2,ai+1,aj+1);
311
      
312
    omdr2obm2  = one - dr2/bm2;
313
    
314
    vbond      = -half*kb*bm2*log(omdr2obm2);
315
    fbond      = -kb/omdr2obm2;
316

    
317
    vtot      += vbond;       /* 35 */
318
    
319
    if (g) {
320
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
321
      ki=IVEC2IS(dt);
322
    }
323
    for (m=0; (m<DIM); m++) {                          /*  15          */
324
      fij=fbond*dx[m];
325
      f[ai][m]+=fij;
326
      f[aj][m]-=fij;
327
      fshift[ki][m]+=fij;
328
      fshift[CENTRAL][m]-=fij;
329
    }
330
  }                                           /*  58 TOTAL    */
331
  return vtot;
332
}
333

    
334
real harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
335
              real *V,real *F)
336
{
337
  const real half=0.5;
338
  real  L1,kk,x0,dx,dx2;
339
  real  v,f,dvdl;
340
  
341
  L1    = 1.0-lambda;
342
  kk    = L1*kA+lambda*kB;
343
  x0    = L1*xA+lambda*xB;
344
  
345
  dx    = x-x0;
346
  dx2   = dx*dx;
347
  
348
  f     = -kk*dx;
349
  v     = half*kk*dx2;
350
  dvdl  = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
351
  
352
  *F    = f;
353
  *V    = v;
354
  
355
  return dvdl;
356
  
357
  /* That was 19 flops */
358
}
359

    
360

    
361
real bonds(int nbonds,
362
           const t_iatom forceatoms[],const t_iparams forceparams[],
363
           const rvec x[],rvec f[],rvec fshift[],
364
           const t_pbc *pbc,const t_graph *g,
365
           real lambda,real *dvdlambda,
366
           const t_mdatoms *md,t_fcdata *fcd)
367
{
368
  int  i,m,ki,ai,aj,type;
369
  real dr,dr2,fbond,vbond,fij,vtot;
370
  rvec dx;
371
  ivec dt;
372

    
373
  vtot = 0.0;
374
  for(i=0; (i<nbonds); ) {
375
    type = forceatoms[i++];
376
    ai   = forceatoms[i++];
377
    aj   = forceatoms[i++];
378
  
379
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);        /*   3                 */
380
    dr2  = iprod(dx,dx);                        /*   5                */
381
    dr   = dr2*invsqrt(dr2);                        /*  10                */
382

    
383
    *dvdlambda += harmonic(forceparams[type].harmonic.krA,
384
                           forceparams[type].harmonic.krB,
385
                           forceparams[type].harmonic.rA,
386
                           forceparams[type].harmonic.rB,
387
                           dr,lambda,&vbond,&fbond);  /*  19  */
388

    
389
    if (dr2 == 0.0)
390
      continue;
391

    
392
    
393
    vtot  += vbond;/* 1*/
394
    fbond *= invsqrt(dr2);                        /*   6                */
395
#ifdef DEBUG
396
    if (debug)
397
      fprintf(debug,"BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
398
              dr,vbond,fbond);
399
#endif
400
    if (g) {
401
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
402
      ki=IVEC2IS(dt);
403
    }
404
    for (m=0; (m<DIM); m++) {                        /*  15                */
405
      fij=fbond*dx[m];
406
      f[ai][m]+=fij;
407
      f[aj][m]-=fij;
408
      fshift[ki][m]+=fij;
409
      fshift[CENTRAL][m]-=fij;
410
    }
411
  }                                        /* 59 TOTAL        */
412
  return vtot;
413
}
414

    
415
real polarize(int nbonds,
416
              const t_iatom forceatoms[],const t_iparams forceparams[],
417
              const rvec x[],rvec f[],rvec fshift[],
418
              const t_pbc *pbc,const t_graph *g,
419
              real lambda,real *dvdlambda,
420
              const t_mdatoms *md,t_fcdata *fcd)
421
{
422
  int  i,m,ki,ai,aj,type;
423
  real dr,dr2,fbond,vbond,fij,vtot,ksh;
424
  rvec dx;
425
  ivec dt;
426

    
427
  vtot = 0.0;
428
  for(i=0; (i<nbonds); ) {
429
    type = forceatoms[i++];
430
    ai   = forceatoms[i++];
431
    aj   = forceatoms[i++];
432
    ksh  = sqr(md->chargeA[aj])*ONE_4PI_EPS0/forceparams[type].polarize.alpha;
433
    if (debug)
434
      fprintf(debug,"POL: ai = %d aj = %d ksh = %.3f\n",ai,aj,ksh);
435
  
436
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);        /*   3                 */
437
    dr2  = iprod(dx,dx);                        /*   5                */
438
    dr   = dr2*invsqrt(dr2);                        /*  10                */
439

    
440
    *dvdlambda += harmonic(ksh,ksh,0,0,dr,lambda,&vbond,&fbond);  /*  19  */
441

    
442
    if (dr2 == 0.0)
443
      continue;
444
    
445
    vtot  += vbond;/* 1*/
446
    fbond *= invsqrt(dr2);                        /*   6                */
447

    
448
    if (g) {
449
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
450
      ki=IVEC2IS(dt);
451
    }
452
    for (m=0; (m<DIM); m++) {                        /*  15                */
453
      fij=fbond*dx[m];
454
      f[ai][m]+=fij;
455
      f[aj][m]-=fij;
456
      fshift[ki][m]+=fij;
457
      fshift[CENTRAL][m]-=fij;
458
    }
459
  }                                        /* 59 TOTAL        */
460
  return vtot;
461
}
462

    
463
real water_pol(int nbonds,
464
               const t_iatom forceatoms[],const t_iparams forceparams[],
465
               const rvec x[],rvec f[],rvec fshift[],
466
               const t_pbc *pbc,const t_graph *g,
467
               real lambda,real *dvdlambda,
468
               const t_mdatoms *md,t_fcdata *fcd)
469
{
470
  /* This routine implements anisotropic polarizibility for water, through
471
   * a shell connected to a dummy with spring constant that differ in the
472
   * three spatial dimensions in the molecular frame.
473
   */
474
  int  i,m,aO,aH1,aH2,aD,aS,type,type0;
475
  rvec dOH1,dOH2,dHH,dOD,dDS,nW,kk,dx,kdx,proj;
476
#ifdef DEBUG
477
  rvec df;
478
#endif
479
  real vtot,fij,r_HH,r_OD,r_nW,tx,ty,tz,qS;
480

    
481
  vtot = 0.0;
482
  if (nbonds > 0) {
483
    type0  = forceatoms[0];
484
    aS     = forceatoms[5];
485
    qS     = md->chargeA[aS];
486
    kk[XX] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_x;
487
    kk[YY] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_y;
488
    kk[ZZ] = sqr(qS)*ONE_4PI_EPS0/forceparams[type0].wpol.al_z;
489
    r_HH   = 1.0/forceparams[type0].wpol.rHH;
490
    r_OD   = 1.0/forceparams[type0].wpol.rOD;
491
    if (debug) {
492
      fprintf(debug,"WPOL: qS  = %10.5f aS = %5d\n",qS,aS);
493
      fprintf(debug,"WPOL: kk  = %10.3f        %10.3f        %10.3f\n",
494
              kk[XX],kk[YY],kk[ZZ]);
495
      fprintf(debug,"WPOL: rOH = %10.3f  rHH = %10.3f  rOD = %10.3f\n",
496
              forceparams[type0].wpol.rOH,
497
              forceparams[type0].wpol.rHH,
498
              forceparams[type0].wpol.rOD);
499
    }
500
    for(i=0; (i<nbonds); i+=6) {
501
      type = forceatoms[i];
502
      if (type != type0)
503
        gmx_fatal(FARGS,"Sorry, type = %d, type0 = %d, file = %s, line = %d",
504
                    type,type0,__FILE__,__LINE__);
505
      aO   = forceatoms[i+1];
506
      aH1  = forceatoms[i+2];
507
      aH2  = forceatoms[i+3];
508
      aD   = forceatoms[i+4];
509
      aS   = forceatoms[i+5];
510
      
511
      /* Compute vectors describing the water frame */
512
      rvec_sub(x[aH1],x[aO], dOH1);
513
      rvec_sub(x[aH2],x[aO], dOH2);
514
      rvec_sub(x[aH2],x[aH1],dHH);
515
      rvec_sub(x[aD], x[aO], dOD);
516
      rvec_sub(x[aS], x[aD], dDS);
517
      oprod(dOH1,dOH2,nW);
518
      
519
      /* Compute inverse length of normal vector 
520
       * (this one could be precomputed, but I'm too lazy now)
521
       */
522
      r_nW = invsqrt(iprod(nW,nW));
523
      /* This is for precision, but does not make a big difference,
524
       * it can go later.
525
       */
526
      r_OD = invsqrt(iprod(dOD,dOD)); 
527
      
528
      /* Normalize the vectors in the water frame */
529
      svmul(r_nW,nW,nW);
530
      svmul(r_HH,dHH,dHH);
531
      svmul(r_OD,dOD,dOD);
532
      
533
      /* Compute displacement of shell along components of the vector */
534
      dx[ZZ] = iprod(dDS,dOD);
535
      /* Compute projection on the XY plane: dDS - dx[ZZ]*dOD */
536
      for(m=0; (m<DIM); m++)
537
        proj[m] = dDS[m]-dx[ZZ]*dOD[m];
538
      
539
      /*dx[XX] = iprod(dDS,nW);
540
        dx[YY] = iprod(dDS,dHH);*/
541
      dx[XX] = iprod(proj,nW);
542
      for(m=0; (m<DIM); m++)
543
        proj[m] -= dx[XX]*nW[m];
544
      dx[YY] = iprod(proj,dHH);
545
      /*#define DEBUG*/
546
#ifdef DEBUG
547
      if (debug) {
548
        fprintf(debug,"WPOL: dx2=%10g  dy2=%10g  dz2=%10g  sum=%10g  dDS^2=%10g\n",
549
                sqr(dx[XX]),sqr(dx[YY]),sqr(dx[ZZ]),iprod(dx,dx),iprod(dDS,dDS));
550
        fprintf(debug,"WPOL: dHH=(%10g,%10g,%10g)\n",dHH[XX],dHH[YY],dHH[ZZ]);
551
        fprintf(debug,"WPOL: dOD=(%10g,%10g,%10g), 1/r_OD = %10g\n",
552
                dOD[XX],dOD[YY],dOD[ZZ],1/r_OD);
553
        fprintf(debug,"WPOL: nW =(%10g,%10g,%10g), 1/r_nW = %10g\n",
554
                nW[XX],nW[YY],nW[ZZ],1/r_nW);
555
        fprintf(debug,"WPOL: dx  =%10g, dy  =%10g, dz  =%10g\n",
556
                dx[XX],dx[YY],dx[ZZ]);
557
         fprintf(debug,"WPOL: dDSx=%10g, dDSy=%10g, dDSz=%10g\n",
558
                dDS[XX],dDS[YY],dDS[ZZ]);
559
      }
560
#endif
561
      /* Now compute the forces and energy */
562
      kdx[XX] = kk[XX]*dx[XX];
563
      kdx[YY] = kk[YY]*dx[YY];
564
      kdx[ZZ] = kk[ZZ]*dx[ZZ];
565
      vtot   += iprod(dx,kdx);
566
      for(m=0; (m<DIM); m++) {
567
        /* This is a tensor operation but written out for speed */
568
        tx        =  nW[m]*kdx[XX];
569
        ty        = dHH[m]*kdx[YY];
570
        tz        = dOD[m]*kdx[ZZ];
571
        fij       = -tx-ty-tz;
572
#ifdef DEBUG
573
        df[m] = fij;
574
#endif
575
        f[aS][m] += fij;
576
        f[aD][m] -= fij;
577
      }
578
#ifdef DEBUG
579
      if (debug) {
580
        fprintf(debug,"WPOL: vwpol=%g\n",0.5*iprod(dx,kdx));
581
        fprintf(debug,"WPOL: df = (%10g, %10g, %10g)\n",df[XX],df[YY],df[ZZ]);
582
      }
583
#endif
584
    }        
585
  }
586
  return 0.5*vtot;
587
}
588

    
589
static real do_1_thole(const rvec xi,const rvec xj,rvec fi,rvec fj,
590
                       const t_pbc *pbc,real qq,
591
                       rvec fshift[],real afac)
592
{
593
  rvec r12;
594
  real r12sq,r12_1,r12n,r12bar,v0,v1,fscal,ebar,fff;
595
  int  m,t;
596
    
597
  t      = pbc_rvec_sub(pbc,xi,xj,r12); /*  3 */
598
  
599
  r12sq  = iprod(r12,r12);              /*  5 */
600
  r12_1  = invsqrt(r12sq);              /*  5 */
601
  r12bar = afac/r12_1;                  /*  5 */
602
  v0     = qq*ONE_4PI_EPS0*r12_1;       /*  2 */
603
  ebar   = exp(-r12bar);                /*  5 */
604
  v1     = (1-(1+0.5*r12bar)*ebar);     /*  4 */
605
  fscal  = ((v0*r12_1)*v1 - v0*afac*ebar*(0.5*r12bar+0.5))*r12_1; /* 9 */
606
  if (debug)
607
    fprintf(debug,"THOLE: v0 = %.3f v1 = %.3f r12= % .3f r12bar = %.3f fscal = %.3f\n",v0,v1,1/r12_1,r12bar,fscal);
608
    
609
  for(m=0; (m<DIM); m++) {
610
    fff    = fscal*r12[m];
611
    fi[m] += fff;
612
    fj[m] -= fff;             
613
    fshift[t][m]       += fff;
614
    fshift[CENTRAL][m] -= fff;
615
  } /* 15 */
616
  
617
  return v0*v1; /* 1 */
618
  /* 54 */
619
}
620

    
621
real thole_pol(int nbonds,
622
               const t_iatom forceatoms[],const t_iparams forceparams[],
623
               const rvec x[],rvec f[],rvec fshift[],
624
               const t_pbc *pbc,const t_graph *g,
625
               real lambda,real *dvdlambda,
626
               const t_mdatoms *md,t_fcdata *fcd)
627
{
628
  /* Interaction between two pairs of particles with opposite charge */
629
  int i,type,a1,da1,a2,da2;
630
  real q1,q2,qq,a,al1,al2,afac;
631
  real V=0;
632
  
633
  for(i=0; (i<nbonds); ) {
634
    type  = forceatoms[i++];
635
    a1    = forceatoms[i++];
636
    da1   = forceatoms[i++];
637
    a2    = forceatoms[i++];
638
    da2   = forceatoms[i++];
639
    q1    = md->chargeA[da1];
640
    q2    = md->chargeA[da2];
641
    a     = forceparams[type].thole.a;
642
    al1   = forceparams[type].thole.alpha1;
643
    al2   = forceparams[type].thole.alpha2;
644
    qq    = q1*q2;
645
    afac  = a*pow(al1*al2,-1.0/6.0);
646
    V += do_1_thole(x[a1], x[a2], f[a1], f[a2], pbc, qq,fshift,afac);
647
    V += do_1_thole(x[da1],x[a2], f[da1],f[a2], pbc,-qq,fshift,afac);
648
    V += do_1_thole(x[a1], x[da2],f[a1], f[da2],pbc,-qq,fshift,afac);
649
    V += do_1_thole(x[da1],x[da2],f[da1],f[da2],pbc, qq,fshift,afac);
650
  }
651
  /* 290 flops */
652
  return V;
653
}
654

    
655
real bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
656
                rvec r_ij,rvec r_kj,real *costh,
657
                int *t1,int *t2)
658
/* Return value is the angle between the bonds i-j and j-k */
659
{
660
  /* 41 FLOPS */
661
  real th;
662
  
663
  *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                        /*  3                */
664
  *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                        /*  3                */
665

    
666
  *costh=cos_angle_non_linear(r_ij,r_kj);                /* 25                */
667

    
668
  th=acos(*costh);                        /* 10                */
669
                                        /* 41 TOTAL        */
670
  return th;
671
}
672

    
673
real angles(int nbonds,
674
            const t_iatom forceatoms[],const t_iparams forceparams[],
675
            const rvec x[],rvec f[],rvec fshift[],
676
            const t_pbc *pbc,const t_graph *g,
677
            real lambda,real *dvdlambda,
678
            const t_mdatoms *md,t_fcdata *fcd)
679
{
680
  int  i,ai,aj,ak,t1,t2,type;
681
  rvec r_ij,r_kj;
682
  real cos_theta,theta,dVdt,va,vtot;
683
  ivec jt,dt_ij,dt_kj;
684
  
685
  vtot = 0.0;
686
  for(i=0; (i<nbonds); ) {
687
    type = forceatoms[i++];
688
    ai   = forceatoms[i++];
689
    aj   = forceatoms[i++];
690
    ak   = forceatoms[i++];
691
    
692
    theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
693
                        r_ij,r_kj,&cos_theta,&t1,&t2);        /*  41                */
694
  
695
    *dvdlambda += harmonic(forceparams[type].harmonic.krA,
696
                           forceparams[type].harmonic.krB,
697
                           forceparams[type].harmonic.rA*DEG2RAD,
698
                           forceparams[type].harmonic.rB*DEG2RAD,
699
                           theta,lambda,&va,&dVdt);  /*  21  */
700
    vtot += va;
701
    
702
    {
703
      int  m;
704
      real snt,st,sth;
705
      real cik,cii,ckk;
706
      real nrkj2,nrij2;
707
      rvec f_i,f_j,f_k;
708
      
709
      snt=sin(theta);                                /*  10                */
710
      if (fabs(snt) < 1e-12)
711
        snt=1e-12;
712
      st  = dVdt/snt;                                /*  10                */
713
      sth = st*cos_theta;                        /*   1                */
714
#ifdef DEBUG
715
      if (debug)
716
        fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
717
                theta*RAD2DEG,va,dVdt);
718
#endif
719
      nrkj2=iprod(r_kj,r_kj);                        /*   5                */
720
      nrij2=iprod(r_ij,r_ij);
721
      
722
      cik=st*invsqrt(nrkj2*nrij2);                /*  12                */ 
723
      cii=sth/nrij2;                                /*  10                */
724
      ckk=sth/nrkj2;                                /*  10                */
725
      
726
      for (m=0; (m<DIM); m++) {                        /*  39                */
727
        f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
728
        f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
729
        f_j[m]=-f_i[m]-f_k[m];
730
        f[ai][m]+=f_i[m];
731
        f[aj][m]+=f_j[m];
732
        f[ak][m]+=f_k[m];
733
      }
734
      if (g) {
735
        copy_ivec(SHIFT_IVEC(g,aj),jt);
736
      
737
        ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
738
        ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
739
        t1=IVEC2IS(dt_ij);
740
        t2=IVEC2IS(dt_kj);
741
      }
742
      rvec_inc(fshift[t1],f_i);
743
      rvec_inc(fshift[CENTRAL],f_j);
744
      rvec_inc(fshift[t2],f_k);
745
    }                                           /* 168 TOTAL        */
746
  }
747
  return vtot;
748
}
749

    
750
real urey_bradley(int nbonds,
751
                  const t_iatom forceatoms[],const t_iparams forceparams[],
752
                  const rvec x[],rvec f[],rvec fshift[],
753
                  const t_pbc *pbc,const t_graph *g,
754
                  real lambda,real *dvdlambda,
755
                  const t_mdatoms *md,t_fcdata *fcd)
756
{
757
  int  i,m,ai,aj,ak,t1,t2,type,ki;
758
  rvec r_ij,r_kj,r_ik;
759
  real cos_theta,theta,dVdt,va,vtot,kth,th0,kUB,r13,dr,dr2,vbond,fbond,fik;
760
  ivec jt,dt_ij,dt_kj,dt_ik;
761
  
762
  vtot = 0.0;
763
  for(i=0; (i<nbonds); ) {
764
    type = forceatoms[i++];
765
    ai   = forceatoms[i++];
766
    aj   = forceatoms[i++];
767
    ak   = forceatoms[i++];
768
    th0  = forceparams[type].u_b.theta*DEG2RAD;
769
    kth  = forceparams[type].u_b.ktheta;
770
    r13  = forceparams[type].u_b.r13;
771
    kUB  = forceparams[type].u_b.kUB;
772
    
773
    theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
774
                        r_ij,r_kj,&cos_theta,&t1,&t2);        /*  41                */
775
  
776
    *dvdlambda += harmonic(kth,kth,th0,th0,theta,lambda,&va,&dVdt);  /*  21  */
777
    vtot += va;
778
    
779
    ki   = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);        /*   3                 */
780
    dr2  = iprod(r_ik,r_ik);                        /*   5                */
781
    dr   = dr2*invsqrt(dr2);                        /*  10                */
782

    
783
    *dvdlambda += harmonic(kUB,kUB,r13,r13,dr,lambda,&vbond,&fbond); /*  19  */
784

    
785
    {
786
      real snt,st,sth;
787
      real cik,cii,ckk;
788
      real nrkj2,nrij2;
789
      rvec f_i,f_j,f_k;
790
      
791
      snt=sin(theta);                                /*  10                */
792
      if (fabs(snt) < 1e-12)
793
        snt=1e-12;
794
      st  = dVdt/snt;                                /*  10                */
795
      sth = st*cos_theta;                        /*   1                */
796
#ifdef DEBUG
797
      if (debug)
798
        fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
799
                theta*RAD2DEG,va,dVdt);
800
#endif
801
      nrkj2=iprod(r_kj,r_kj);                        /*   5                */
802
      nrij2=iprod(r_ij,r_ij);
803
      
804
      cik=st*invsqrt(nrkj2*nrij2);                /*  12                */ 
805
      cii=sth/nrij2;                                /*  10                */
806
      ckk=sth/nrkj2;                                /*  10                */
807
      
808
      for (m=0; (m<DIM); m++) {                        /*  39                */
809
        f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
810
        f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
811
        f_j[m]=-f_i[m]-f_k[m];
812
        f[ai][m]+=f_i[m];
813
        f[aj][m]+=f_j[m];
814
        f[ak][m]+=f_k[m];
815
      }
816
      if (g) {
817
        copy_ivec(SHIFT_IVEC(g,aj),jt);
818
      
819
        ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
820
        ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
821
        t1=IVEC2IS(dt_ij);
822
        t2=IVEC2IS(dt_kj);
823
      }
824
      rvec_inc(fshift[t1],f_i);
825
      rvec_inc(fshift[CENTRAL],f_j);
826
      rvec_inc(fshift[t2],f_k);
827
    }                                           /* 168 TOTAL        */
828
    /* Time for the bond calculations */
829
    if (dr2 == 0.0)
830
      continue;
831

    
832
    vtot  += vbond;  /* 1*/
833
    fbond *= invsqrt(dr2);                        /*   6                */
834
    
835
    if (g) {
836
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,ak),dt_ik);
837
      ki=IVEC2IS(dt_ik);
838
    }
839
    for (m=0; (m<DIM); m++) {                        /*  15                */
840
      fik=fbond*r_ik[m];
841
      f[ai][m]+=fik;
842
      f[ak][m]-=fik;
843
      fshift[ki][m]+=fik;
844
      fshift[CENTRAL][m]-=fik;
845
    }
846
  }
847
  return vtot;
848
}
849

    
850
real quartic_angles(int nbonds,
851
                    const t_iatom forceatoms[],const t_iparams forceparams[],
852
                    const rvec x[],rvec f[],rvec fshift[],
853
                    const t_pbc *pbc,const t_graph *g,
854
                    real lambda,real *dvdlambda,
855
                    const t_mdatoms *md,t_fcdata *fcd)
856
{
857
  int  i,j,ai,aj,ak,t1,t2,type;
858
  rvec r_ij,r_kj;
859
  real cos_theta,theta,dt,dVdt,va,dtp,c,vtot;
860
  ivec jt,dt_ij,dt_kj;
861
  
862
  vtot = 0.0;
863
  for(i=0; (i<nbonds); ) {
864
    type = forceatoms[i++];
865
    ai   = forceatoms[i++];
866
    aj   = forceatoms[i++];
867
    ak   = forceatoms[i++];
868

    
869
    theta  = bond_angle(x[ai],x[aj],x[ak],pbc,
870
                        r_ij,r_kj,&cos_theta,&t1,&t2);        /*  41                */
871

    
872
    dt = theta - forceparams[type].qangle.theta*DEG2RAD; /* 2          */
873

    
874
    dVdt = 0;
875
    va = forceparams[type].qangle.c[0];
876
    dtp = 1.0;
877
    for(j=1; j<=4; j++) {
878
      c = forceparams[type].qangle.c[j];
879
      dVdt -= j*c*dtp;
880
      dtp *= dt;
881
      va += c*dtp;
882
    }
883
    /* 20 */
884

    
885
    vtot += va;
886
    
887
    {
888
      int  m;
889
      real snt,st,sth;
890
      real cik,cii,ckk;
891
      real nrkj2,nrij2;
892
      rvec f_i,f_j,f_k;
893
      
894
      snt=sin(theta);                                /*  10                */
895
      if (fabs(snt) < 1e-12)
896
        snt=1e-12;
897
      st  = dVdt/snt;                                /*  10                */
898
      sth = st*cos_theta;                        /*   1                */
899
#ifdef DEBUG
900
      if (debug)
901
        fprintf(debug,"ANGLES: theta = %10g  vth = %10g  dV/dtheta = %10g\n",
902
                theta*RAD2DEG,va,dVdt);
903
#endif
904
      nrkj2=iprod(r_kj,r_kj);                        /*   5                */
905
      nrij2=iprod(r_ij,r_ij);
906
      
907
      cik=st*invsqrt(nrkj2*nrij2);                /*  12                */ 
908
      cii=sth/nrij2;                                /*  10                */
909
      ckk=sth/nrkj2;                                /*  10                */
910
      
911
      for (m=0; (m<DIM); m++) {                        /*  39                */
912
        f_i[m]=-(cik*r_kj[m]-cii*r_ij[m]);
913
        f_k[m]=-(cik*r_ij[m]-ckk*r_kj[m]);
914
        f_j[m]=-f_i[m]-f_k[m];
915
        f[ai][m]+=f_i[m];
916
        f[aj][m]+=f_j[m];
917
        f[ak][m]+=f_k[m];
918
      }
919
      if (g) {
920
        copy_ivec(SHIFT_IVEC(g,aj),jt);
921
      
922
        ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
923
        ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
924
        t1=IVEC2IS(dt_ij);
925
        t2=IVEC2IS(dt_kj);
926
      }
927
      rvec_inc(fshift[t1],f_i);
928
      rvec_inc(fshift[CENTRAL],f_j);
929
      rvec_inc(fshift[t2],f_k);
930
    }                                           /* 160 TOTAL        */
931
  }
932
  return vtot;
933
}
934

    
935
real dih_angle(const rvec xi,const rvec xj,const rvec xk,const rvec xl,
936
               const t_pbc *pbc,
937
               rvec r_ij,rvec r_kj,rvec r_kl,rvec m,rvec n,
938
               real *cos_phi,real *sign,int *t1,int *t2,int *t3)
939
{
940
  real ipr,phi;
941

    
942
  *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                       /*  3                 */
943
  *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                        /*  3                */
944
  *t3 = pbc_rvec_sub(pbc,xk,xl,r_kl);                        /*  3                */
945

    
946
  oprod(r_ij,r_kj,m);                         /*  9                 */
947
  oprod(r_kj,r_kl,n);                        /*  9                */
948
  *cos_phi=cos_angle_non_linear(m,n);                 /* 41                 */
949
  phi=acos(*cos_phi);                         /* 10                 */
950
  ipr=iprod(r_ij,n);                         /*  5                 */
951
  (*sign)=(ipr<0.0)?-1.0:1.0;
952
  phi=(*sign)*phi;                         /*  1                */
953
                                        /* 84 TOTAL        */
954
  return phi;
955
}
956

    
957
void do_dih_fup(int i,int j,int k,int l,real ddphi,
958
                rvec r_ij,rvec r_kj,rvec r_kl,
959
                rvec m,rvec n,rvec f[],rvec fshift[],
960
                const t_pbc *pbc,const t_graph *g,
961
                const rvec x[],int t1,int t2,int t3)
962
{
963
  /* 143 FLOPS */
964
  rvec f_i,f_j,f_k,f_l;
965
  rvec uvec,vvec,svec,dx_jl;
966
  real iprm,iprn,nrkj,nrkj2;
967
  real a,p,q;
968
  ivec jt,dt_ij,dt_kj,dt_lj;  
969
  
970
  iprm  = iprod(m,m);                /*  5         */
971
  iprn  = iprod(n,n);                /*  5        */
972
  if ((iprm > GMX_REAL_EPS) && (iprn > GMX_REAL_EPS)) {
973
    nrkj2 = iprod(r_kj,r_kj);        /*  5        */
974
    nrkj  = nrkj2*invsqrt(nrkj2);        /* 10        */
975
    a     = -ddphi*nrkj/iprm;        /* 11        */
976
    svmul(a,m,f_i);                /*  3        */
977
    a     = ddphi*nrkj/iprn;        /* 11        */
978
    svmul(a,n,f_l);                /*  3         */
979
    p     = iprod(r_ij,r_kj);        /*  5        */
980
    p    /= nrkj2;                /* 10        */
981
    q     = iprod(r_kl,r_kj);        /*  5        */
982
    q    /= nrkj2;                /* 10        */
983
    svmul(p,f_i,uvec);                /*  3        */
984
    svmul(q,f_l,vvec);                /*  3        */
985
    rvec_sub(uvec,vvec,svec);        /*  3        */
986
    rvec_sub(f_i,svec,f_j);        /*  3        */
987
    rvec_add(f_l,svec,f_k);        /*  3        */
988
    rvec_inc(f[i],f_i);           /*  3        */
989
    rvec_dec(f[j],f_j);           /*  3        */
990
    rvec_dec(f[k],f_k);           /*  3        */
991
    rvec_inc(f[l],f_l);           /*  3        */
992
    
993
    if (g) {
994
      copy_ivec(SHIFT_IVEC(g,j),jt);
995
      ivec_sub(SHIFT_IVEC(g,i),jt,dt_ij);
996
      ivec_sub(SHIFT_IVEC(g,k),jt,dt_kj);
997
      ivec_sub(SHIFT_IVEC(g,l),jt,dt_lj);
998
      t1=IVEC2IS(dt_ij);
999
      t2=IVEC2IS(dt_kj);
1000
      t3=IVEC2IS(dt_lj);
1001
    }    
1002
    else
1003
      t3 = pbc_rvec_sub(pbc,x[l],x[j],dx_jl);
1004
    
1005
    rvec_inc(fshift[t1],f_i);
1006
    rvec_dec(fshift[CENTRAL],f_j);
1007
    rvec_dec(fshift[t2],f_k);
1008
    rvec_inc(fshift[t3],f_l);
1009
  }
1010
  /* 112 TOTAL         */
1011
}
1012

    
1013

    
1014
real dopdihs(real cpA,real cpB,real phiA,real phiB,int mult,
1015
             real phi,real lambda,real *V,real *F)
1016
{
1017
  real v,dvdl,mdphi,v1,sdphi,ddphi;
1018
  real L1   = 1.0 - lambda;
1019
  real ph0  = (L1*phiA + lambda*phiB)*DEG2RAD;
1020
  real dph0 = (phiB - phiA)*DEG2RAD;
1021
  real cp   = L1*cpA + lambda*cpB;
1022
  
1023
  mdphi =  mult*phi - ph0;
1024
  sdphi = sin(mdphi);
1025
  ddphi = -cp*mult*sdphi;
1026
  v1    = 1.0 + cos(mdphi);
1027
  v     = cp*v1;
1028
  
1029
  dvdl  = (cpB - cpA)*v1 + cp*dph0*sdphi;
1030
  
1031
  *V = v;
1032
  *F = ddphi;
1033
  
1034
  return dvdl;
1035
  
1036
  /* That was 40 flops */
1037
}
1038

    
1039
static real dopdihs_min(real cpA,real cpB,real phiA,real phiB,int mult,
1040
                       real phi,real lambda,real *V,real *F)
1041
     /* similar to dopdihs, except for a minus sign  *
1042
      * and a different treatment of mult/phi0       */
1043
{
1044
  real v,dvdl,mdphi,v1,sdphi,ddphi;
1045
  real L1   = 1.0 - lambda;
1046
  real ph0  = (L1*phiA + lambda*phiB)*DEG2RAD;
1047
  real dph0 = (phiB - phiA)*DEG2RAD;
1048
  real cp   = L1*cpA + lambda*cpB;
1049
  
1050
  mdphi = mult*(phi-ph0);
1051
  sdphi = sin(mdphi);
1052
  ddphi = cp*mult*sdphi;
1053
  v1    = 1.0-cos(mdphi);
1054
  v     = cp*v1;
1055
  
1056
  dvdl  = (cpB-cpA)*v1 + cp*dph0*sdphi;
1057
  
1058
  *V = v;
1059
  *F = ddphi;
1060
  
1061
  return dvdl;
1062
  
1063
  /* That was 40 flops */
1064
}
1065

    
1066
real pdihs(int nbonds,
1067
           const t_iatom forceatoms[],const t_iparams forceparams[],
1068
           const rvec x[],rvec f[],rvec fshift[],
1069
           const t_pbc *pbc,const t_graph *g,
1070
           real lambda,real *dvdlambda,
1071
           const t_mdatoms *md,t_fcdata *fcd)
1072
{
1073
  int  i,type,ai,aj,ak,al;
1074
  int  t1,t2,t3;
1075
  rvec r_ij,r_kj,r_kl,m,n;
1076
  real phi,cos_phi,sign,ddphi,vpd,vtot;
1077

    
1078
  vtot = 0.0;
1079
  for(i=0; (i<nbonds); ) {
1080
    type = forceatoms[i++];
1081
    ai   = forceatoms[i++];
1082
    aj   = forceatoms[i++];
1083
    ak   = forceatoms[i++];
1084
    al   = forceatoms[i++];
1085
    
1086
    phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1087
                  &cos_phi,&sign,&t1,&t2,&t3);                        /*  84                 */
1088
                
1089
    *dvdlambda += dopdihs(forceparams[type].pdihs.cpA,
1090
                          forceparams[type].pdihs.cpB,
1091
                          forceparams[type].pdihs.phiA,
1092
                          forceparams[type].pdihs.phiB,
1093
                          forceparams[type].pdihs.mult,
1094
                          phi,lambda,&vpd,&ddphi);
1095
                       
1096
    vtot += vpd;
1097
    do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1098
               f,fshift,pbc,g,x,t1,t2,t3);                        /* 112                */
1099

    
1100
#ifdef DEBUG
1101
    fprintf(debug,"pdih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
1102
            ai,aj,ak,al,cos_phi,phi);
1103
#endif
1104
  } /* 223 TOTAL         */
1105

    
1106
  return vtot;
1107
}
1108

    
1109

    
1110

    
1111
real idihs(int nbonds,
1112
           const t_iatom forceatoms[],const t_iparams forceparams[],
1113
           const rvec x[],rvec f[],rvec fshift[],
1114
           const t_pbc *pbc,const t_graph *g,
1115
           real lambda,real *dvdlambda,
1116
           const t_mdatoms *md,t_fcdata *fcd)
1117
{
1118
  int  i,type,ai,aj,ak,al;
1119
  int  t1,t2,t3;
1120
  real phi,phi0,dphi0,cos_phi,ddphi,sign,vtot;
1121
  rvec r_ij,r_kj,r_kl,m,n;
1122
  real L1,kk,dp,dp2,kA,kB,pA,pB,dvdl;
1123

    
1124
  L1 = 1.0-lambda;
1125
  dvdl = 0;
1126

    
1127
  vtot = 0.0;
1128
  for(i=0; (i<nbonds); ) {
1129
    type = forceatoms[i++];
1130
    ai   = forceatoms[i++];
1131
    aj   = forceatoms[i++];
1132
    ak   = forceatoms[i++];
1133
    al   = forceatoms[i++];
1134
    
1135
    phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1136
                  &cos_phi,&sign,&t1,&t2,&t3);                        /*  84                */
1137
    
1138
    /* phi can jump if phi0 is close to Pi/-Pi, which will cause huge
1139
     * force changes if we just apply a normal harmonic.
1140
     * Instead, we first calculate phi-phi0 and take it modulo (-Pi,Pi).
1141
     * This means we will never have the periodicity problem, unless
1142
     * the dihedral is Pi away from phiO, which is very unlikely due to
1143
     * the potential.
1144
     */
1145
    kA = forceparams[type].harmonic.krA;
1146
    kB = forceparams[type].harmonic.krB;
1147
    pA = forceparams[type].harmonic.rA;
1148
    pB = forceparams[type].harmonic.rB;
1149

    
1150
    kk    = L1*kA + lambda*kB;
1151
    phi0  = (L1*pA + lambda*pB)*DEG2RAD;
1152
    dphi0 = (pB - pA)*DEG2RAD;
1153

    
1154
    /* dp = (phi-phi0), modulo (-pi,pi) */
1155
    dp = phi-phi0;  
1156
    /* dp cannot be outside (-2*pi,2*pi) */
1157
    if (dp >= M_PI)
1158
      dp -= 2*M_PI;
1159
    else if(dp < -M_PI)
1160
      dp += 2*M_PI;
1161
    
1162
    dp2 = dp*dp;
1163

    
1164
    vtot += 0.5*kk*dp2;
1165
    ddphi = -kk*dp;
1166
    
1167
    dvdl += 0.5*(kB - kA)*dp2 - kk*dphi0*dp;
1168

    
1169
    do_dih_fup(ai,aj,ak,al,(real)(-ddphi),r_ij,r_kj,r_kl,m,n,
1170
               f,fshift,pbc,g,x,t1,t2,t3);                        /* 112                */
1171
    /* 217 TOTAL        */
1172
#ifdef DEBUG
1173
    if (debug)
1174
      fprintf(debug,"idih: (%d,%d,%d,%d) cp=%g, phi=%g\n",
1175
              ai,aj,ak,al,cos_phi,phi);
1176
#endif
1177
  }
1178
  
1179
  *dvdlambda += dvdl;
1180
  return vtot;
1181
}
1182

    
1183

    
1184
real posres(int nbonds,
1185
            const t_iatom forceatoms[],const t_iparams forceparams[],
1186
            const rvec x[],rvec f[],rvec fshift[],
1187
            const t_pbc *pbc,const t_graph *g,
1188
            real lambda,real *dvdlambda,
1189
            const t_mdatoms *md,t_fcdata *fcd)
1190
{
1191
  int  i,ai,m,type,ki;
1192
  const t_iparams *pr;
1193
  real v,vtot,fm,*fc;
1194
  rvec dx;
1195

    
1196
  vtot = 0.0;
1197
  for(i=0; (i<nbonds); ) {
1198
    type = forceatoms[i++];
1199
    ai   = forceatoms[i++];
1200
    pr   = &forceparams[type];
1201
    
1202
    /*
1203
    if (pbc == NULL)
1204
      rvec_sub(x[ai],forceparams[type].posres.pos0A,dx);
1205
    else
1206
    pbc_dx(pbc,x[ai],forceparams[type].posres.pos0A,dx);*/
1207
    
1208
    ki = pbc_rvec_sub(pbc,x[ai],forceparams[type].posres.pos0A,dx);
1209
    
1210
    v=0;
1211
    for (m=0; (m<DIM); m++) {
1212
      *dvdlambda += harmonic(pr->posres.fcA[m],pr->posres.fcB[m],
1213
                             0,pr->posres.pos0B[m]-pr->posres.pos0A[m],
1214
                             dx[m],lambda,&v,&fm);
1215
      vtot += v;
1216
      f[ai][m] += fm;
1217
      /* Position restraints should not contribute to the virial.
1218
       * We need to subtract this contribution !!!
1219
       * Needs to be implemented, but requires an extra virial matrix.
1220
       */
1221
    }
1222
  }
1223

    
1224
  return vtot;
1225
}
1226

    
1227
static real low_angres(int nbonds,
1228
                       const t_iatom forceatoms[],const t_iparams forceparams[],
1229
                       const rvec x[],rvec f[],rvec fshift[],
1230
                       const t_pbc *pbc,const t_graph *g,
1231
                       real lambda,real *dvdlambda,
1232
                       bool bZAxis)
1233
{
1234
  int  i,m,type,ai,aj,ak,al;
1235
  int  t1,t2;
1236
  real phi,cos_phi,vid,vtot,dVdphi;
1237
  rvec r_ij,r_kl,f_i,f_k;
1238
  real st,sth,sin_phi,nrij2,nrkl2,c,cij,ckl;
1239

    
1240
  ivec dt;  
1241
  t2 = 0; /* avoid warning with gcc-3.3. It is never used uninitialized */
1242

    
1243
  vtot = 0.0;
1244
  ak=al=0; /* to avoid warnings */
1245
  for(i=0; i<nbonds; ) {
1246
    type = forceatoms[i++];
1247
    ai   = forceatoms[i++];
1248
    aj   = forceatoms[i++];
1249
    t1   = pbc_rvec_sub(pbc,x[aj],x[ai],r_ij);                    /*  3                */
1250
    if (!bZAxis) {      
1251
      ak   = forceatoms[i++];
1252
      al   = forceatoms[i++];
1253
      t2   = pbc_rvec_sub(pbc,x[al],x[ak],r_kl);           /*  3                */
1254
    } else {
1255
      r_kl[XX] = 0;
1256
      r_kl[YY] = 0;
1257
      r_kl[ZZ] = 1;
1258
    }
1259

    
1260
    cos_phi = cos_angle_non_linear(r_ij,r_kl);                /* 25                */
1261
    phi     = acos(cos_phi);                    /* 10           */
1262

    
1263
    *dvdlambda += dopdihs_min(forceparams[type].pdihs.cpA,
1264
                              forceparams[type].pdihs.cpB,
1265
                              forceparams[type].pdihs.phiA,
1266
                              forceparams[type].pdihs.phiB,
1267
                              forceparams[type].pdihs.mult,
1268
                              phi,lambda,&vid,&dVdphi); /*  40  */
1269
    
1270
    vtot += vid;
1271

    
1272
    sin_phi = sin(phi);                                /*  10                */
1273
    if (fabs(sin_phi) < 1e-12)
1274
      sin_phi=1e-12;
1275
    st  = -dVdphi/sin_phi;                        /*  10                */
1276
    sth = st*cos_phi;                                /*   1                */
1277
    nrij2 = iprod(r_ij,r_ij);                        /*   5                */
1278
    nrkl2 = iprod(r_kl,r_kl);                   /*   5          */
1279
    
1280
    c   = st*invsqrt(nrij2*nrkl2);                /*  11                */ 
1281
    cij = sth/nrij2;                                /*  10                */
1282
    ckl = sth/nrkl2;                                /*  10                */
1283
    
1284
    for (m=0; m<DIM; m++) {                        /*  18+18       */
1285
      f_i[m] = (c*r_kl[m]-cij*r_ij[m]);
1286
      f[ai][m] += f_i[m];
1287
      f[aj][m] -= f_i[m];
1288
      if (!bZAxis) {
1289
        f_k[m] = (c*r_ij[m]-ckl*r_kl[m]);
1290
        f[ak][m] += f_k[m];
1291
        f[al][m] -= f_k[m];
1292
      }
1293
    }
1294
    
1295
    if (g) {
1296
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1297
      t1=IVEC2IS(dt);
1298
    }
1299
    rvec_inc(fshift[t1],f_i);
1300
    rvec_dec(fshift[CENTRAL],f_i);
1301
    if (!bZAxis) {
1302
      if (g) {
1303
        ivec_sub(SHIFT_IVEC(g,ak),SHIFT_IVEC(g,al),dt);
1304
        t2=IVEC2IS(dt);
1305
      }
1306
      rvec_inc(fshift[t2],f_k);
1307
      rvec_dec(fshift[CENTRAL],f_k);
1308
    }
1309
  }
1310

    
1311
  return vtot;  /*  191 / 164 (bZAxis)  total  */
1312
}
1313

    
1314
real angres(int nbonds,
1315
            const t_iatom forceatoms[],const t_iparams forceparams[],
1316
            const rvec x[],rvec f[],rvec fshift[],
1317
            const t_pbc *pbc,const t_graph *g,
1318
            real lambda,real *dvdlambda,
1319
            const t_mdatoms *md,t_fcdata *fcd)
1320
{
1321
  return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1322
                    lambda,dvdlambda,FALSE);
1323
}
1324

    
1325
real angresz(int nbonds,
1326
             const t_iatom forceatoms[],const t_iparams forceparams[],
1327
             const rvec x[],rvec f[],rvec fshift[],
1328
             const t_pbc *pbc,const t_graph *g,
1329
             real lambda,real *dvdlambda,
1330
             const t_mdatoms *md,t_fcdata *fcd)
1331
{
1332
  return low_angres(nbonds,forceatoms,forceparams,x,f,fshift,pbc,g,
1333
                    lambda,dvdlambda,TRUE);
1334
}
1335

    
1336

    
1337
real unimplemented(int nbonds,
1338
                   const t_iatom forceatoms[],const t_iparams forceparams[],
1339
                   const rvec x[],rvec f[],rvec fshift[],
1340
                   const t_pbc *pbc,const t_graph *g,
1341
                   real lambda,real *dvdlambda,
1342
                   const t_mdatoms *md,t_fcdata *fcd)
1343
{
1344
  gmx_impl("*** you are using a not implemented function");
1345

    
1346
  return 0.0; /* To make the compiler happy */
1347
}
1348

    
1349
real rbdihs(int nbonds,
1350
            const t_iatom forceatoms[],const t_iparams forceparams[],
1351
            const rvec x[],rvec f[],rvec fshift[],
1352
            const t_pbc *pbc,const t_graph *g,
1353
            real lambda,real *dvdlambda,
1354
            const t_mdatoms *md,t_fcdata *fcd)
1355
{
1356
  const real c0=0.0,c1=1.0,c2=2.0,c3=3.0,c4=4.0,c5=5.0;
1357
  int  type,ai,aj,ak,al,i,j;
1358
  int  t1,t2,t3;
1359
  rvec r_ij,r_kj,r_kl,m,n;
1360
  real parmA[NR_RBDIHS];
1361
  real parmB[NR_RBDIHS];
1362
  real parm[NR_RBDIHS];
1363
  real phi,cos_phi,rbp,rbpBA;
1364
  real v,sign,ddphi,sin_phi;
1365
  real cosfac,vtot;
1366
  real L1   = 1.0-lambda;
1367
  real dvdl=0;
1368

    
1369
  vtot = 0.0;
1370
  for(i=0; (i<nbonds); ) {
1371
    type = forceatoms[i++];
1372
    ai   = forceatoms[i++];
1373
    aj   = forceatoms[i++];
1374
    ak   = forceatoms[i++];
1375
    al   = forceatoms[i++];
1376

    
1377
    phi=dih_angle(x[ai],x[aj],x[ak],x[al],pbc,r_ij,r_kj,r_kl,m,n,
1378
                  &cos_phi,&sign,&t1,&t2,&t3);                        /*  84                */
1379

    
1380
    /* Change to polymer convention */
1381
    if (phi < c0)
1382
      phi += M_PI;
1383
    else
1384
      phi -= M_PI;                        /*   1                */
1385
    cos_phi = -cos_phi;                                        /*   1                */
1386
    
1387
    sin_phi=sin(phi);
1388

    
1389
    for(j=0; (j<NR_RBDIHS); j++) {
1390
      parmA[j] = forceparams[type].rbdihs.rbcA[j];
1391
      parmB[j] = forceparams[type].rbdihs.rbcB[j];
1392
      parm[j]  = L1*parmA[j]+lambda*parmB[j];
1393
    }
1394
    /* Calculate cosine powers */
1395
    /* Calculate the energy */
1396
    /* Calculate the derivative */
1397

    
1398
    v       = parm[0];
1399
    dvdl   += (parmB[0]-parmA[0]);
1400
    ddphi   = c0;
1401
    cosfac  = c1;
1402
    
1403
    rbp     = parm[1];
1404
    rbpBA   = parmB[1]-parmA[1];
1405
    ddphi  += rbp*cosfac;
1406
    cosfac *= cos_phi;
1407
    v      += cosfac*rbp;
1408
    dvdl   += cosfac*rbpBA;
1409
    rbp     = parm[2];
1410
    rbpBA   = parmB[2]-parmA[2];    
1411
    ddphi  += c2*rbp*cosfac;
1412
    cosfac *= cos_phi;
1413
    v      += cosfac*rbp;
1414
    dvdl   += cosfac*rbpBA;
1415
    rbp     = parm[3];
1416
    rbpBA   = parmB[3]-parmA[3];
1417
    ddphi  += c3*rbp*cosfac;
1418
    cosfac *= cos_phi;
1419
    v      += cosfac*rbp;
1420
    dvdl   += cosfac*rbpBA;
1421
    rbp     = parm[4];
1422
    rbpBA   = parmB[4]-parmA[4];
1423
    ddphi  += c4*rbp*cosfac;
1424
    cosfac *= cos_phi;
1425
    v      += cosfac*rbp;
1426
    dvdl   += cosfac*rbpBA;
1427
    rbp     = parm[5];
1428
    rbpBA   = parmB[5]-parmA[5];
1429
    ddphi  += c5*rbp*cosfac;
1430
    cosfac *= cos_phi;
1431
    v      += cosfac*rbp;
1432
    dvdl   += cosfac*rbpBA;
1433
   
1434
    ddphi = -ddphi*sin_phi;                                /*  11                */
1435
    
1436
    do_dih_fup(ai,aj,ak,al,ddphi,r_ij,r_kj,r_kl,m,n,
1437
               f,fshift,pbc,g,x,t1,t2,t3);                /* 112                */
1438
    vtot += v;
1439
  }  
1440
  *dvdlambda += dvdl;
1441

    
1442
  return vtot;
1443
}
1444

    
1445
/***********************************************************
1446
 *
1447
 *   G R O M O S  9 6   F U N C T I O N S
1448
 *
1449
 ***********************************************************/
1450
real g96harmonic(real kA,real kB,real xA,real xB,real x,real lambda,
1451
                 real *V,real *F)
1452
{
1453
  const real half=0.5;
1454
  real  L1,kk,x0,dx,dx2;
1455
  real  v,f,dvdl;
1456
  
1457
  L1    = 1.0-lambda;
1458
  kk    = L1*kA+lambda*kB;
1459
  x0    = L1*xA+lambda*xB;
1460
  
1461
  dx    = x-x0;
1462
  dx2   = dx*dx;
1463
  
1464
  f     = -kk*dx;
1465
  v     = half*kk*dx2;
1466
  dvdl  = half*(kB-kA)*dx2 + (xA-xB)*kk*dx;
1467
  
1468
  *F    = f;
1469
  *V    = v;
1470
  
1471
  return dvdl;
1472
  
1473
  /* That was 21 flops */
1474
}
1475

    
1476
real g96bonds(int nbonds,
1477
              const t_iatom forceatoms[],const t_iparams forceparams[],
1478
              const rvec x[],rvec f[],rvec fshift[],
1479
              const t_pbc *pbc,const t_graph *g,
1480
              real lambda,real *dvdlambda,
1481
              const t_mdatoms *md,t_fcdata *fcd)
1482
{
1483
  int  i,m,ki,ai,aj,type;
1484
  real dr2,fbond,vbond,fij,vtot;
1485
  rvec dx;
1486
  ivec dt;
1487
  
1488
  vtot = 0.0;
1489
  for(i=0; (i<nbonds); ) {
1490
    type = forceatoms[i++];
1491
    ai   = forceatoms[i++];
1492
    aj   = forceatoms[i++];
1493
  
1494
    ki   = pbc_rvec_sub(pbc,x[ai],x[aj],dx);                /*   3                 */
1495
    dr2  = iprod(dx,dx);                                /*   5                */
1496
      
1497
    *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
1498
                              forceparams[type].harmonic.krB,
1499
                              forceparams[type].harmonic.rA,
1500
                              forceparams[type].harmonic.rB,
1501
                              dr2,lambda,&vbond,&fbond);
1502

    
1503
    vtot  += 0.5*vbond;                             /* 1*/
1504
#ifdef DEBUG
1505
    if (debug)
1506
      fprintf(debug,"G96-BONDS: dr = %10g  vbond = %10g  fbond = %10g\n",
1507
              sqrt(dr2),vbond,fbond);
1508
#endif
1509
   
1510
    if (g) {
1511
      ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);
1512
      ki=IVEC2IS(dt);
1513
    }
1514
    for (m=0; (m<DIM); m++) {                        /*  15                */
1515
      fij=fbond*dx[m];
1516
      f[ai][m]+=fij;
1517
      f[aj][m]-=fij;
1518
      fshift[ki][m]+=fij;
1519
      fshift[CENTRAL][m]-=fij;
1520
    }
1521
  }                                        /* 44 TOTAL        */
1522
  return vtot;
1523
}
1524

    
1525
real g96bond_angle(const rvec xi,const rvec xj,const rvec xk,const t_pbc *pbc,
1526
                   rvec r_ij,rvec r_kj,
1527
                   int *t1,int *t2)
1528
/* Return value is the angle between the bonds i-j and j-k */
1529
{
1530
  real costh;
1531
  
1532
  *t1 = pbc_rvec_sub(pbc,xi,xj,r_ij);                        /*  3                */
1533
  *t2 = pbc_rvec_sub(pbc,xk,xj,r_kj);                        /*  3                */
1534

    
1535
  costh=cos_angle_non_linear(r_ij,r_kj);                /* 25                */
1536
                                        /* 41 TOTAL        */
1537
  return costh;
1538
}
1539

    
1540
real g96angles(int nbonds,
1541
               const t_iatom forceatoms[],const t_iparams forceparams[],
1542
               const rvec x[],rvec f[],rvec fshift[],
1543
               const t_pbc *pbc,const t_graph *g,
1544
               real lambda,real *dvdlambda,
1545
               const t_mdatoms *md,t_fcdata *fcd)
1546
{
1547
  int  i,ai,aj,ak,type,m,t1,t2;
1548
  rvec r_ij,r_kj;
1549
  real cos_theta,dVdt,va,vtot;
1550
  real rij_1,rij_2,rkj_1,rkj_2,rijrkj_1;
1551
  rvec f_i,f_j,f_k;
1552
  ivec jt,dt_ij,dt_kj;
1553
  
1554
  vtot = 0.0;
1555
  for(i=0; (i<nbonds); ) {
1556
    type = forceatoms[i++];
1557
    ai   = forceatoms[i++];
1558
    aj   = forceatoms[i++];
1559
    ak   = forceatoms[i++];
1560
    
1561
    cos_theta  = g96bond_angle(x[ai],x[aj],x[ak],pbc,r_ij,r_kj,&t1,&t2);
1562

    
1563
    *dvdlambda += g96harmonic(forceparams[type].harmonic.krA,
1564
                              forceparams[type].harmonic.krB,
1565
                              forceparams[type].harmonic.rA,
1566
                              forceparams[type].harmonic.rB,
1567
                              cos_theta,lambda,&va,&dVdt);
1568
    vtot    += va;
1569
    
1570
    rij_1    = invsqrt(iprod(r_ij,r_ij));
1571
    rkj_1    = invsqrt(iprod(r_kj,r_kj));
1572
    rij_2    = rij_1*rij_1;
1573
    rkj_2    = rkj_1*rkj_1;
1574
    rijrkj_1 = rij_1*rkj_1;                     /* 23 */
1575
    
1576
#ifdef DEBUG
1577
    if (debug)
1578
      fprintf(debug,"G96ANGLES: costheta = %10g  vth = %10g  dV/dct = %10g\n",
1579
              cos_theta,va,dVdt);
1580
#endif
1581
    for (m=0; (m<DIM); m++) {                        /*  42        */
1582
      f_i[m]=dVdt*(r_kj[m]*rijrkj_1 - r_ij[m]*rij_2*cos_theta);
1583
      f_k[m]=dVdt*(r_ij[m]*rijrkj_1 - r_kj[m]*rkj_2*cos_theta);
1584
      f_j[m]=-f_i[m]-f_k[m];
1585
      f[ai][m]+=f_i[m];
1586
      f[aj][m]+=f_j[m];
1587
      f[ak][m]+=f_k[m];
1588
    }
1589
    
1590
    if (g) {
1591
      copy_ivec(SHIFT_IVEC(g,aj),jt);
1592
      
1593
      ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
1594
      ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
1595
      t1=IVEC2IS(dt_ij);
1596
      t2=IVEC2IS(dt_kj);
1597
    }      
1598
    rvec_inc(fshift[t1],f_i);
1599
    rvec_inc(fshift[CENTRAL],f_j);
1600
    rvec_inc(fshift[t2],f_k);               /* 9 */
1601
    /* 163 TOTAL        */
1602
  }
1603
  return vtot;
1604
}
1605

    
1606
real cross_bond_bond(int nbonds,
1607
                     const t_iatom forceatoms[],const t_iparams forceparams[],
1608
                     const rvec x[],rvec f[],rvec fshift[],
1609
                     const t_pbc *pbc,const t_graph *g,
1610
                     real lambda,real *dvdlambda,
1611
                     const t_mdatoms *md,t_fcdata *fcd)
1612
{
1613
  /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
1614
   * pp. 842-847
1615
   */
1616
  int  i,ai,aj,ak,type,m,t1,t2;
1617
  rvec r_ij,r_kj;
1618
  real vtot,vrr,s1,s2,r1,r2,r1e,r2e,krr;
1619
  rvec f_i,f_j,f_k;
1620
  ivec jt,dt_ij,dt_kj;
1621
  
1622
  vtot = 0.0;
1623
  for(i=0; (i<nbonds); ) {
1624
    type = forceatoms[i++];
1625
    ai   = forceatoms[i++];
1626
    aj   = forceatoms[i++];
1627
    ak   = forceatoms[i++];
1628
    r1e  = forceparams[type].cross_bb.r1e;
1629
    r2e  = forceparams[type].cross_bb.r2e;
1630
    krr  = forceparams[type].cross_bb.krr;
1631
    
1632
    /* Compute distance vectors ... */
1633
    t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
1634
    t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
1635
    
1636
    /* ... and their lengths */
1637
    r1 = norm(r_ij);
1638
    r2 = norm(r_kj);
1639
    
1640
    /* Deviations from ideality */
1641
    s1 = r1-r1e;
1642
    s2 = r2-r2e;
1643
    
1644
    /* Energy (can be negative!) */
1645
    vrr   = krr*s1*s2;
1646
    vtot += vrr;
1647
    
1648
    /* Forces */
1649
    svmul(-krr*s2/r1,r_ij,f_i);
1650
    svmul(-krr*s1/r2,r_kj,f_k);
1651
    
1652
    for (m=0; (m<DIM); m++) {                        /*  12        */
1653
      f_j[m]    = -f_i[m] - f_k[m];
1654
      f[ai][m] += f_i[m];
1655
      f[aj][m] += f_j[m];
1656
      f[ak][m] += f_k[m];
1657
    }
1658
    
1659
    /* Virial stuff */
1660
    if (g) {
1661
      copy_ivec(SHIFT_IVEC(g,aj),jt);
1662
      
1663
      ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
1664
      ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
1665
      t1=IVEC2IS(dt_ij);
1666
      t2=IVEC2IS(dt_kj);
1667
    }      
1668
    rvec_inc(fshift[t1],f_i);
1669
    rvec_inc(fshift[CENTRAL],f_j);
1670
    rvec_inc(fshift[t2],f_k);               /* 9 */
1671
    /* 163 TOTAL        */
1672
  }
1673
  return vtot;
1674
}
1675

    
1676
real cross_bond_angle(int nbonds,
1677
                      const t_iatom forceatoms[],const t_iparams forceparams[],
1678
                      const rvec x[],rvec f[],rvec fshift[],
1679
                      const t_pbc *pbc,const t_graph *g,
1680
                      real lambda,real *dvdlambda,
1681
                      const t_mdatoms *md,t_fcdata *fcd)
1682
{
1683
  /* Potential from Lawrence and Skimmer, Chem. Phys. Lett. 372 (2003)
1684
   * pp. 842-847
1685
   */
1686
  int  i,ai,aj,ak,type,m,t1,t2,t3;
1687
  rvec r_ij,r_kj,r_ik;
1688
  real vtot,vrt,s1,s2,s3,r1,r2,r3,r1e,r2e,r3e,krt,k1,k2,k3;
1689
  rvec f_i,f_j,f_k;
1690
  ivec jt,dt_ij,dt_kj;
1691
  
1692
  vtot = 0.0;
1693
  for(i=0; (i<nbonds); ) {
1694
    type = forceatoms[i++];
1695
    ai   = forceatoms[i++];
1696
    aj   = forceatoms[i++];
1697
    ak   = forceatoms[i++];
1698
    r1e  = forceparams[type].cross_ba.r1e;
1699
    r2e  = forceparams[type].cross_ba.r2e;
1700
    r3e  = forceparams[type].cross_ba.r3e;
1701
    krt  = forceparams[type].cross_ba.krt;
1702
    
1703
    /* Compute distance vectors ... */
1704
    t1 = pbc_rvec_sub(pbc,x[ai],x[aj],r_ij);
1705
    t2 = pbc_rvec_sub(pbc,x[ak],x[aj],r_kj);
1706
    t3 = pbc_rvec_sub(pbc,x[ai],x[ak],r_ik);
1707
    
1708
    /* ... and their lengths */
1709
    r1 = norm(r_ij);
1710
    r2 = norm(r_kj);
1711
    r3 = norm(r_ik);
1712
    
1713
    /* Deviations from ideality */
1714
    s1 = r1-r1e;
1715
    s2 = r2-r2e;
1716
    s3 = r3-r3e;
1717
    
1718
    /* Energy (can be negative!) */
1719
    vrt   = krt*s3*(s1+s2);
1720
    vtot += vrt;
1721
    
1722
    /* Forces */
1723
    k1 = -krt*(s3/r1);
1724
    k2 = -krt*(s3/r2);
1725
    k3 = -krt*(s1+s2)/r3;
1726
    for(m=0; (m<DIM); m++) {
1727
      f_i[m] = k1*r_ij[m] + k3*r_ik[m];
1728
      f_k[m] = k2*r_kj[m] - k3*r_ik[m];
1729
      f_j[m] = -f_i[m] - f_k[m];
1730
    }
1731
    
1732
    for (m=0; (m<DIM); m++) {                        /*  12        */
1733
      f[ai][m] += f_i[m];
1734
      f[aj][m] += f_j[m];
1735
      f[ak][m] += f_k[m];
1736
    }
1737
    
1738
    /* Virial stuff */
1739
    if (g) {
1740
      copy_ivec(SHIFT_IVEC(g,aj),jt);
1741
      
1742
      ivec_sub(SHIFT_IVEC(g,ai),jt,dt_ij);
1743
      ivec_sub(SHIFT_IVEC(g,ak),jt,dt_kj);
1744
      t1=IVEC2IS(dt_ij);
1745
      t2=IVEC2IS(dt_kj);
1746
    }      
1747
    rvec_inc(fshift[t1],f_i);
1748
    rvec_inc(fshift[CENTRAL],f_j);
1749
    rvec_inc(fshift[t2],f_k);               /* 9 */
1750
    /* 163 TOTAL        */
1751
  }
1752
  return vtot;
1753
}
1754