src/mdlib/fshake.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 fshake(iatom,ncon,nit,maxnit,
  48      $     dist2,xp,rij,m2,omega,invmass,tt,lagr,error)
  49
  50       implicit none
  51
  52       integer*4 iatom(*)
  53       integer*4 ncon,nit,maxnit
  54       integer*4 error
  55       real*4    dist2(*),xp(*),rij(*),m2(*),invmass(*),tt(*),lagr(*)
  56       real*4    omega
  57
  58       integer*4 ll,i,j,i3,j3,l3,nconv,iconv
  59       integer*4 ix,iy,iz,jx,jy,jz
  60       real*4  toler,rpij2,rrpr,tx,ty,tz,diff,acor,im,jm
  61       real*4  xh,yh,zh,rijx,rijy,rijz
  62       real*4  mytol
  63
  64       parameter(mytol=1e-6)
  65
  66       error=0
  67       do nit=1,maxnit
  68          nconv=0
  69          do ll=1,ncon
  70             l3    = 3*(ll-1)
  71             rijx  = rij(l3+1)
  72             rijy  = rij(l3+2)
  73             rijz  = rij(l3+3)
  74             i     = iatom(l3+2)
  75             j     = iatom(l3+3)
  76             i3    = 3*i
  77             j3    = 3*j
  78             ix    = i3+1
  79             iy    = i3+2
  80             iz    = i3+3
  81             jx    = j3+1
  82             jy    = j3+2
  83             jz    = j3+3
  84
  85             tx      = xp(ix)-xp(jx)
  86             ty      = xp(iy)-xp(jy)
  87             tz      = xp(iz)-xp(jz)
  88             rpij2   = tx*tx+ty*ty+tz*tz
  89
  90             toler   = dist2(ll)
  91             diff    = toler-rpij2
  92             iconv   = abs(diff)*tt(ll)
  93
  94             if (iconv .ne. 0) then
  95                nconv   = nconv + iconv
  96                rrpr    = rijx*tx+rijy*ty+rijz*tz
  97
  98                if (rrpr .lt. mytol*toler) then
  99                   error   = ll
 100                else
 101                   acor    = omega*diff*m2(ll)/rrpr
 102                   lagr(ll) = lagr(ll) + acor
 103                   xh      = rijx*acor
 104                   yh      = rijy*acor
 105                   zh      = rijz*acor
 106                   im      = invmass(i+1)
 107                   jm      = invmass(j+1)
 108                   if ((im .ne. 0) .and. (jm .ne. 0)) then
 109                      xp(ix) = xp(ix) + xh*im
 110                      xp(iy) = xp(iy) + yh*im
 111                      xp(iz) = xp(iz) + zh*im
 112                      xp(jx) = xp(jx) - xh*jm
 113                      xp(jy) = xp(jy) - yh*jm
 114                      xp(jz) = xp(jz) - zh*jm
 115                   else
 116                      if ((im .eq. 0) .and. (jm .ne. 0)) then
 117                         xp(ix) = xp(ix) + xh*jm
 118                         xp(iy) = xp(iy) + yh*jm
 119                         xp(iz) = xp(iz) + zh*jm
 120                      else if ((jm .eq. 0) .and. (im .ne. 0)) then
 121                         xp(jx) = xp(jx) - xh*im
 122                         xp(jy) = xp(jy) - yh*im
 123                         xp(jz) = xp(jz) - zh*im
 124                      else
 125                         error = ll
 126                      endif
 127                   endif
 128                end if
 129             end if
 130          end do
 131
 132          if (nconv .eq. 0) goto 10
 133          if (error .ne. 0) goto 10
 134       end do
 135
 136  10   return
 137
 138       end
 139