src/mdlib/fsettle.f

   1
   2 C
   3 C                This source code is part of
   4 C
   5 C                 G   R   O   M   A   C   S
   6 C
   7 C          GROningen MAchine for Chemical Simulations
   8 C
   9 C                        VERSION 3.0
  10 C
  11 C Copyright (c) 1991-2001
  12 C BIOSON Research Institute, Dept. of Biophysical Chemistry
  13 C University of Groningen, The Netherlands
  14 C
  15 C This program is free software; you can redistribute it and/or
  16 C modify it under the terms of the GNU General Public License
  17 C as published by the Free Software Foundation; either version 2
  18 C of the License, or (at your option) any later version.
  19 C
  20 C If you want to redistribute modifications, please consider that
  21 C scientific software is very special. Version control is crucial -
  22 C bugs must be traceable. We will be happy to consider code for
  23 C inclusion in the official distribution, but derived work must not
  24 C be called official GROMACS. Details are found in the README & COPYING
  25 C files - if they are missing, get the official version at www.gromacs.org.
  26 C
  27 C To help us fund GROMACS development, we humbly ask that you cite
  28 C the papers on the package - you can find them in the top README file.
  29 C
  30 C Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  31 C
  32 C And Hey:
  33 C GROup of MAchos and Cynical Suckers
  34
  35 c     IMPORTANT IMPORTANT IMPORTANT IMPORTANT!
  36 c     Note that this file comes in two flavours -
  37 c     fshake.f for single precision and fshaked.f
  38 c     for double precision. The only difference is
  39 c     the size of the real variables.
  40 c     This is an unfortunate, but necessary setup
  41 c     since not all f77 compilers (e.g. g77) have
  42 c     switches to change the real size, and neither
  43 c     do all f77 compilers support preprocessing.
  44 c     Thus, if you edit one of the files, make sure
  45 c      to change to other similarly!
  46
  47       subroutine fsettle(nshake,owptr,b4,after,
  48      &     dOH,dHH,mO,mH,error)
  49       implicit none
  50 c*****************************************************************
  51 c                                                               **
  52 c    Subroutine : setlep - reset positions of TIP3P waters      **
  53 c    Author : Shuichi Miyamoto                                  **
  54 c    Date of last update : Oct. 1, 1992                         **
  55 c                                                               **
  56 c    Reference for the SETTLE algorithm                         **
  57 c           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). **
  58 c                                                               **
  59 c*****************************************************************
  60       integer nshake,owptr(*),error
  61
  62       real*4  b4(*),after(*),mO,mH
  63      $     ,dOH,dHH
  64       real*4  wo,wh,wohh,ra,rb,rc,rc2,tmp,tmp2
  65       real*4 gama, beta, alpa, xcom, ycom, zcom, al2be2
  66       real*4 axlng, aylng, azlng, trns11, trns21, trns31, trns12,
  67      &     trns22, trns32, trns13, trns23, trns33, cosphi,
  68      &     costhe, sinphi, sinthe, cospsi, xaksxd, yaksxd,
  69      &     xakszd, yakszd, zakszd, zaksxd, xaksyd,
  70      &     xb0, yb0, zb0, xc0, yc0, zc0, xa1
  71       real*4 ya1, za1, xb1, yb1
  72       real*4 zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi,
  73      &     xa3, ya3, za3, xb3, yb3, zb3, xc3, yc3, zc3,
  74      &     xb0d, yb0d, xc0d, yc0d, xa1d, ya1d,
  75      &     za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d,
  76      &     xb2d, yb2d, yc2d,
  77      &     xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d
  78
  79       integer i,ow1,hw2,hw3
  80
  81       error= -1
  82       wo   = mO
  83       wh   = mH
  84       wohh = mO+2.0*mH
  85       rc   = dHH/2.0
  86       ra   = 2.0*wh*sqrt(dOH*dOH-rc*rc)/wohh
  87       rb   = sqrt(dOH*dOH-rc*rc)-ra
  88       rc2  = dHH
  89       wo   = wo/wohh
  90       wh   = wh/wohh
  91 c
  92 cray compiler directive ignore vector dependencies
  93 c$dir ivdep
  94       do i=1,nshake
  95 c                                                --- Step1  A1' ---
  96 c
  97          ow1 = 3*owptr(i)+1
  98          hw2 = ow1+3
  99          hw3 = ow1+6
 100          xb0 = b4(hw2)   - b4(ow1)
 101          yb0 = b4(hw2+1) - b4(ow1+1)
 102          zb0 = b4(hw2+2) - b4(ow1+2)
 103          xc0 = b4(hw3)   - b4(ow1)
 104          yc0 = b4(hw3+1) - b4(ow1+1)
 105          zc0 = b4(hw3+2) - b4(ow1+2)
 106 c
 107          xcom = ( after(ow1  )*wo +
 108      &        (after(hw2  ) + after(hw3  )) * wh )
 109          ycom = ( after(ow1+1)*wo +
 110      &        (after(hw2+1) + after(hw3+1)) * wh )
 111          zcom = ( after(ow1+2)*wo +
 112      &        (after(hw2+2) + after(hw3+2)) * wh )
 113 c
 114          xa1 = after(ow1  ) - xcom
 115          ya1 = after(ow1+1) - ycom
 116          za1 = after(ow1+2) - zcom
 117          xb1 = after(hw2  ) - xcom
 118          yb1 = after(hw2+1) - ycom
 119          zb1 = after(hw2+2) - zcom
 120          xc1 = after(hw3  ) - xcom
 121          yc1 = after(hw3+1) - ycom
 122          zc1 = after(hw3+2) - zcom
 123 c
 124          xaksZd = yb0*zc0 - zb0*yc0
 125          yaksZd = zb0*xc0 - xb0*zc0
 126          zaksZd = xb0*yc0 - yb0*xc0
 127          xaksXd = ya1*zaksZd - za1*yaksZd
 128          yaksXd = za1*xaksZd - xa1*zaksZd
 129          zaksXd = xa1*yaksZd - ya1*xaksZd
 130          xaksYd = yaksZd*zaksXd - zaksZd*yaksXd
 131          yaksYd = zaksZd*xaksXd - xaksZd*zaksXd
 132          zaksYd = xaksZd*yaksXd - yaksZd*xaksXd
 133 c
 134          axlng = 1.0/sqrt ( xaksXd * xaksXd + yaksXd * yaksXd
 135      &                                  + zaksXd * zaksXd )
 136          aylng = 1.0/sqrt ( xaksYd * xaksYd + yaksYd * yaksYd
 137      &                                  + zaksYd * zaksYd )
 138          azlng = 1.0/sqrt ( xaksZd * xaksZd + yaksZd * yaksZd
 139      &                                  + zaksZd * zaksZd )
 140          trns11 = xaksXd * axlng
 141          trns21 = yaksXd * axlng
 142          trns31 = zaksXd * axlng
 143          trns12 = xaksYd * aylng
 144          trns22 = yaksYd * aylng
 145          trns32 = zaksYd * aylng
 146          trns13 = xaksZd * azlng
 147          trns23 = yaksZd * azlng
 148          trns33 = zaksZd * azlng
 149 c
 150          xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0
 151          yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0
 152          xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0
 153          yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0
 154          xa1d = trns11*xa1 + trns21*ya1 + trns31*za1
 155          ya1d = trns12*xa1 + trns22*ya1 + trns32*za1
 156          za1d = trns13*xa1 + trns23*ya1 + trns33*za1
 157          xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1
 158          yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1
 159          zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1
 160          xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1
 161          yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1
 162          zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1
 163 c
 164          sinphi = za1d / ra
 165          tmp    = 1.0 - sinphi*sinphi
 166          if ( tmp .le. 0.0) then
 167             error = i-1
 168             cosphi = 0
 169          else
 170             cosphi = sqrt (tmp)
 171          endif
 172          sinpsi = ( zb1d - zc1d ) / (rc2 * cosphi)
 173          tmp2   = 1.0 - sinpsi*sinpsi
 174          if ( tmp2 .le. 0.0 ) then
 175             error = i-1
 176             cospsi = 0
 177          else
 178             cospsi = sqrt (tmp2)
 179          endif
 180 c
 181          ya2d =   ra * cosphi
 182          xb2d = - rc * cospsi
 183 c        xc2d =   rc * cospsi
 184          yb2d = - rb * cosphi - rc *sinpsi * sinphi
 185          yc2d = - rb * cosphi + rc *sinpsi * sinphi
 186 c        xb2d2 = xb2d * xb2d
 187 c        hh2 = 4.d0 * xb2d2 + (yb2d-yc2d) * (yb2d-yc2d)
 188 c    &                      + (zb1d-zc1d) * (zb1d-zc1d)
 189 c        deltx = 2.d0 * xb2d + sqrt ( 4.d0 * xb2d2 - hh2 + hhhh )
 190 c        xb2d = xb2d - deltx * 0.5d0
 191 c        xc2d = xc2d + deltx * 0.5d0
 192 c
 193 c                                                --- Step3  al,be,ga ---
 194 c
 195          alpa = ( xb2d * (xb0d-xc0d) + yb0d * yb2d + yc0d * yc2d )
 196          beta = ( xb2d * (yc0d-yb0d) + xb0d * yb2d + xc0d * yc2d )
 197          gama = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d
 198 c
 199          al2be2 = alpa * alpa + beta * beta
 200          sinthe = ( alpa*gama - beta * sqrt ( al2be2 - gama * gama ) )
 201      &            / al2be2
 202 c
 203 c                                                --- Step4  A3' ---
 204 c
 205          costhe = sqrt (1.d0 - sinthe * sinthe )
 206          xa3d = - ya2d * sinthe
 207          ya3d =   ya2d * costhe
 208          za3d = za1d
 209          xb3d =   xb2d * costhe - yb2d * sinthe
 210          yb3d =   xb2d * sinthe + yb2d * costhe
 211          zb3d = zb1d
 212          xc3d = - xb2d * costhe - yc2d * sinthe
 213          yc3d = - xb2d * sinthe + yc2d * costhe
 214          zc3d = zc1d
 215 c                                                --- Step5  A3 ---
 216 c
 217 c
 218          xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d
 219          ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d
 220          za3 = trns31*xa3d + trns32*ya3d + trns33*za3d
 221          xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d
 222          yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d
 223          zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d
 224          xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d
 225          yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d
 226          zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d
 227 c
 228          after(ow1  ) = xcom + xa3
 229          after(ow1+1) = ycom + ya3
 230          after(ow1+2) = zcom + za3
 231          after(hw2  ) = xcom + xb3
 232          after(hw2+1) = ycom + yb3
 233          after(hw2+2) = zcom + zb3
 234          after(hw3  ) = xcom + xc3
 235          after(hw3+1) = ycom + yc3
 236          after(hw3+2) = zcom + zc3
 237 c
 238       end do
 239       return
 240       end
 241