src/contrib/intest.f

   1 C
   2 C                This source code is part of
   3 C
   4 C                 G   R   O   M   A   C   S
   5 C
   6 C          GROningen MAchine for Chemical Simulations
   7 C
   8 C                        VERSION 3.0
   9 C
  10 C Copyright (c) 1991-2001
  11 C BIOSON Research Institute, Dept. of Biophysical Chemistry
  12 C University of Groningen, The Netherlands
  13 C
  14 C This program is free software; you can redistribute it and/or
  15 C modify it under the terms of the GNU General Public License
  16 C as published by the Free Software Foundation; either version 2
  17 C of the License, or (at your option) any later version.
  18 C
  19 C If you want to redistribute modifications, please consider that
  20 C scientific software is very special. Version control is crucial -
  21 C bugs must be traceable. We will be happy to consider code for
  22 C inclusion in the official distribution, but derived work must not
  23 C be called official GROMACS. Details are found in the README & COPYING
  24 C files - if they are missing, get the official version at www.gromacs.org.
  25 C
  26 C To help us fund GROMACS development, we humbly ask that you cite
  27 C the papers on the package - you can find them in the top README file.
  28 C
  29 C Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  30 C
  31 C And Hey:
  32 C GROup of MAchos and Cynical Suckers
  33
  34 C
  35 C     This code is meant to be called from C routines.
  36 C     Therefore all indices start at 0, although the arrays
  37 C     start at 1, if an array contains an index we must add 1 to it.
  38 C     EG: jjnr points to particles starting at 0
  39 C         type is indexed from 1 to ...
  40 C
  41
  42       subroutine FORLJC(ix,iy,iz,qi,
  43      $     pos,nj,type,jjnr,charge,nbfp,
  44      $     faction,fip,
  45      $     Vc,Vnb)
  46
  47       implicit none
  48
  49       real      ix,iy,iz,qi
  50       real      pos(*),charge(*),faction(*),fip(3)
  51       integer*4 nj,jjnr(*),type(*)
  52       real      Vc,Vnb,nbfp(*)
  53
  54       integer   k,jnr,j3,tj
  55       real      twelve,six
  56       real      fX,fY,fZ
  57       real      rijX,rijY,rijZ
  58       real      fijscal,vijcoul
  59       real      vctot,vnbtot
  60       real      rinv1,rinv2,rinv6
  61       real      fjx,fjy,fjz
  62       real      tx,ty,tz,vnb6,vnb12
  63
  64       parameter(twelve=12.0,six=6.0)
  65
  66       fX     = 0
  67       fY     = 0
  68       fZ     = 0
  69       vctot  = 0
  70       vnbtot = 0
  71
  72 cray compiler directive ignore vector dependencies
  73 c$dir ivdep
  74       do k=1,nj
  75          jnr   = jjnr(k)+1
  76          j3    = 3*jnr-2
  77          rijX  = ix - pos(j3)
  78          rijY  = iy - pos(j3+1)
  79          rijZ  = iz - pos(j3+2)
  80
  81          rinv1       = 1.0/sqrt((rijX*rijX)+(rijY*rijY)+(rijZ*rijZ))
  82          rinv2       = rinv1*rinv1
  83          rinv6       = rinv2*rinv2*rinv2
  84
  85          tj          = 2*type(jnr)+1
  86          vnb6        = nbfp(tj)*rinv6
  87          vnb12       = nbfp(tj+1)*rinv6*rinv6
  88          vijcoul     = qi*charge(jnr)*rinv1
  89
  90          vctot       = vctot+vijcoul
  91          vnbtot      = vnbtot+vnb12-vnb6
  92          fijscal     = (twelve*vnb12-six*vnb6+vijcoul)*rinv2
  93
  94          fjx           = faction(j3)
  95          tx            = rijX*fijscal
  96          fX            = fX + tx
  97          faction(j3)   = fjx - tx
  98          fjy           = faction(j3+1)
  99          ty            = rijY*fijscal
 100          fY            = fY + ty
 101          faction(j3+1) = fjy - ty
 102          fjz           = faction(j3+2)
 103          tz            = rijZ*fijscal
 104          fZ            = fZ + tz
 105          faction(j3+2) = fjz - tz
 106
 107       end do
 108
 109       fip(1) = fX
 110       fip(2) = fY
 111       fip(3) = fZ
 112       Vc     = Vc  + vctot
 113       Vnb    = Vnb + vnbtot
 114
 115       return
 116
 117       end
 118
 119       program main
 120
 121       implicit none
 122
 123       integer maxatom,maxx,maxlist,maxtype
 124       parameter(maxatom=1000,maxx=3*maxatom,maxlist=100)
 125       parameter(maxtype=19)
 126
 127       real*4 ix,iy,iz,qi,pos(maxx),faction(maxx),fip(3)
 128       real*4 charge(maxatom),nbfp(2*maxtype),Vc,Vnb
 129       integer type(maxatom),jjnr(maxlist)
 130       integer i,i3,j
 131
 132 c     setup benchmark
 133       do i=1,maxtype
 134          nbfp(2*i-1) = 1e-3
 135          nbfp(2*i)   = 1e-6
 136       end do
 137
 138       type(1)=1
 139       do i=2,maxatom
 140          type(i)=1+mod(type(i-1)+91,maxtype)
 141       end do
 142
 143       do i=1,maxatom
 144          i3=3*(i-1)
 145          pos(i3+1) = i
 146          pos(i3+2) = i
 147          pos(i3+3) = 1
 148
 149          charge(i) = mod(i,2)-0.5
 150       end do
 151
 152       jjnr(1) = 13
 153       do i=2,maxlist
 154          jjnr(i) = mod(jjnr(i-1)+13,maxatom)
 155       end do
 156
 157       ix = 0.0
 158       iy = 0.0
 159       iz = 0.0
 160       qi = 1.0
 161       Vc = 0.0
 162       Vnb= 0.0
 163
 164 c     run it
 165
 166
 167       do j=1,100
 168          do i=1,maxatom
 169
 170             call FORLJC(ix,iy,iz,qi,
 171      $           pos,maxlist,type,jjnr,charge,nbfp,
 172      $           faction,fip,
 173      $           Vc,Vnb)
 174
 175          end do
 176       end do
 177
 178       print *,Vc, Vnb
 179
 180       stop
 181       end