src/mdlib/csettle.c

   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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Getting the Right Output Means no Artefacts in Calculating Stuff
  31  */
  32 static char *SRCID_csettle_c = "$Id$";
  33 #include <math.h>
  34 #include <stdio.h>
  35 #include "vec.h"
  36 #include "constr.h"
  37 #include "fatal.h"
  38
  39 #ifdef DEBUG
  40 static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3)
  41 {
  42   rvec dOH1,dOH2,dHH;
  43   int  m;
  44
  45   for(m=0; (m<DIM); m++) {
  46     dOH1[m]=x[OW1+m]-x[HW2+m];
  47     dOH2[m]=x[OW1+m]-x[HW3+m];
  48     dHH[m]=x[HW2+m]-x[HW3+m];
  49   }
  50   fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d,  dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
  51           title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH));
  52 }
  53 #endif
  54
  55 void csettle(FILE *fp,int nshake, int owptr[],real b4[], real after[],
  56              real dOH,real dHH,real mO,real mH,int *error)
  57 {
  58   /* ***************************************************************** */
  59   /*                                                               ** */
  60   /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
  61   /*    Author : Shuichi Miyamoto                                  ** */
  62   /*    Date of last update : Oct. 1, 1992                         ** */
  63   /*                                                               ** */
  64   /*    Reference for the SETTLE algorithm                         ** */
  65   /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
  66   /*                                                               ** */
  67   /* ***************************************************************** */
  68
  69   /* Initialized data */
  70   static bool bFirst=TRUE;
  71   static real wo,wh,wohh;
  72   static real ra,rb,rc,rc2,rone;
  73 #ifdef DEBUG_PRES
  74   static int step = 0;
  75 #endif
  76
  77
  78   /* Local variables */
  79   real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
  80   real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
  81     trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
  82     cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
  83     xb0, yb0, zb0, xc0, yc0, zc0, xa1;
  84   real ya1, za1, xb1, yb1;
  85   real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
  86     xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
  87     za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
  88     xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
  89   real t1,t2;
  90
  91   int i, ow1, hw2, hw3;
  92
  93   *error=-1;
  94   if (bFirst) {
  95     fprintf(fp,"Going to use C-settle (%d waters)\n",nshake);
  96     wo     = mO;
  97     wh     = mH;
  98     wohh   = mO+2.0*mH;
  99     rc     = dHH/2.0;
 100     ra     = 2.0*wh*sqrt(dOH*dOH-rc*rc)/wohh;
 101     rb     = sqrt(dOH*dOH-rc*rc)-ra;
 102     rc2    = dHH;
 103     rone   = 1.0;
 104
 105     wo    /= wohh;
 106     wh    /= wohh;
 107
 108     fprintf(fp,"wo = %g, wh =%g, wohh = %g, rc = %g, ra = %g\n",
 109             wo,wh,wohh,rc,ra);
 110     fprintf(fp,"rb = %g, rc2 = %g, rone = %g, dHH = %g, dOH = %g\n",
 111             rb,rc2,rone,dHH,dOH);
 112
 113     bFirst = FALSE;
 114   }
 115 #ifdef PRAGMAS
 116 #pragma ivdep
 117 #endif
 118   for (i = 0; i < nshake; ++i) {
 119     /*    --- Step1  A1' ---      */
 120     ow1 = owptr[i] * 3;
 121     hw2 = ow1 + 3;
 122     hw3 = ow1 + 6;
 123     xb0 = b4[hw2    ] - b4[ow1];
 124     yb0 = b4[hw2 + 1] - b4[ow1 + 1];
 125     zb0 = b4[hw2 + 2] - b4[ow1 + 2];
 126     xc0 = b4[hw3    ] - b4[ow1];
 127     yc0 = b4[hw3 + 1] - b4[ow1 + 1];
 128     zc0 = b4[hw3 + 2] - b4[ow1 + 2];
 129     /* 6 flops */
 130
 131     xcom = (after[ow1    ] * wo + (after[hw2    ] + after[hw3    ]) * wh);
 132     ycom = (after[ow1 + 1] * wo + (after[hw2 + 1] + after[hw3 + 1]) * wh);
 133     zcom = (after[ow1 + 2] * wo + (after[hw2 + 2] + after[hw3 + 2]) * wh);
 134     /* 12 flops */
 135
 136     xa1 = after[ow1    ] - xcom;
 137     ya1 = after[ow1 + 1] - ycom;
 138     za1 = after[ow1 + 2] - zcom;
 139     xb1 = after[hw2    ] - xcom;
 140     yb1 = after[hw2 + 1] - ycom;
 141     zb1 = after[hw2 + 2] - zcom;
 142     xc1 = after[hw3    ] - xcom;
 143     yc1 = after[hw3 + 1] - ycom;
 144     zc1 = after[hw3 + 2] - zcom;
 145     /* 9 flops */
 146
 147     xakszd = yb0 * zc0 - zb0 * yc0;
 148     yakszd = zb0 * xc0 - xb0 * zc0;
 149     zakszd = xb0 * yc0 - yb0 * xc0;
 150     xaksxd = ya1 * zakszd - za1 * yakszd;
 151     yaksxd = za1 * xakszd - xa1 * zakszd;
 152     zaksxd = xa1 * yakszd - ya1 * xakszd;
 153     xaksyd = yakszd * zaksxd - zakszd * yaksxd;
 154     yaksyd = zakszd * xaksxd - xakszd * zaksxd;
 155     zaksyd = xakszd * yaksxd - yakszd * xaksxd;
 156     /* 27 flops */
 157
 158     axlng = invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
 159     aylng = invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
 160     azlng = invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
 161
 162     trns11 = xaksxd * axlng;
 163     trns21 = yaksxd * axlng;
 164     trns31 = zaksxd * axlng;
 165     trns12 = xaksyd * aylng;
 166     trns22 = yaksyd * aylng;
 167     trns32 = zaksyd * aylng;
 168     trns13 = xakszd * azlng;
 169     trns23 = yakszd * azlng;
 170     trns33 = zakszd * azlng;
 171     /* 24 flops */
 172
 173     xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
 174     yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
 175     xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
 176     yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
 177     /*
 178     xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 179     ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 180     */
 181     za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
 182     xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
 183     yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
 184     zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
 185     xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
 186     yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
 187     zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
 188     /* 65 flops */
 189
 190     sinphi = za1d / ra;
 191     tmp    = rone - sinphi * sinphi;
 192     if (tmp <= 0) {
 193       *error = i;
 194       cosphi = 0;
 195     }
 196     else
 197       cosphi = sqrt(tmp);
 198     sinpsi = (zb1d - zc1d) / (rc2 * cosphi);
 199     tmp2   = rone - sinpsi * sinpsi;
 200     if (tmp2 <= 0) {
 201       *error = i;
 202       cospsi = 0;
 203     }
 204     else
 205       cospsi = sqrt(tmp2);
 206     /* 46 flops */
 207
 208     ya2d =  ra * cosphi;
 209     xb2d = -rc * cospsi;
 210     t1   = -rb * cosphi;
 211     t2   =  rc * sinpsi * sinphi;
 212     yb2d =  t1 - t2;
 213     yc2d =  t1 + t2;
 214     /* 7 flops */
 215
 216     /*     --- Step3  al,be,ga                --- */
 217     alpa   = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
 218     beta   = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
 219     gama   = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
 220     al2be2 = alpa * alpa + beta * beta;
 221     sinthe = (alpa * gama - beta * sqrt(al2be2 - gama * gama)) / al2be2;
 222     /* 47 flops */
 223
 224     /*  --- Step4  A3' --- */
 225     costhe = sqrt(rone - sinthe * sinthe);
 226     xa3d = -ya2d * sinthe;
 227     ya3d = ya2d * costhe;
 228     za3d = za1d;
 229     xb3d = xb2d * costhe - yb2d * sinthe;
 230     yb3d = xb2d * sinthe + yb2d * costhe;
 231     zb3d = zb1d;
 232     xc3d = -xb2d * costhe - yc2d * sinthe;
 233     yc3d = -xb2d * sinthe + yc2d * costhe;
 234     zc3d = zc1d;
 235     /* 26 flops */
 236
 237     /*    --- Step5  A3 --- */
 238     xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
 239     ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
 240     za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
 241     xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
 242     yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
 243     zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
 244     xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
 245     yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
 246     zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
 247     /* 45 flops */
 248
 249     after[ow1] = xcom + xa3;
 250     after[ow1 + 1] = ycom + ya3;
 251     after[ow1 + 2] = zcom + za3;
 252     after[hw2] = xcom + xb3;
 253     after[hw2 + 1] = ycom + yb3;
 254     after[hw2 + 2] = zcom + zb3;
 255     after[hw3] = xcom + xc3;
 256     after[hw3 + 1] = ycom + yc3;
 257     after[hw3 + 2] = zcom + zc3;
 258     /* 9 flops */
 259 #ifdef DEBUG
 260     check_cons(fp,"settle",after,ow1,hw2,hw3);
 261 #endif
 262   }
 263 }
 264