src/contrib/intest.f

   1 C
   2 C        @(#) in_loopf.f 1.17 18 Aug 1996
   3 C
   4 C        This  source-code  is  part  of
   5 C
   6 C        G    R    O    M    A    C    S
   7 C
   8 C  GROningen MAchine for Chemical Simulations
   9 C
  10 C  Copyright (c) 1990-1995,
  11 C  BIOSON Research Institute, Dept. of Biophysical Chemistry,
  12 C  University of Groningen, The Netherlands
  13 C
  14 C  Please refer to:
  15 C  GROMACS: A Message Passing Parallel Molecular Dynamics Implementation
  16 C  H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  17 C  Comp. Phys. Comm. 1995 (in press)
  18 C
  19 C  Also check out our WWW page:
  20 C  http://rugmd0.chem.rug.nl/~gmx/gmx.cgi
  21 C  or e-mail to:
  22 C  gromacs@chem.rug.nl
  23 C
  24 C  And Hey:
  25 C  Gyas ROwers Mature At Cryogenic Speed
  26 C
  27
  28 C
  29 C     This code is meant to be called from C routines.
  30 C     Therefore all indices start at 0, although the arrays
  31 C     start at 1, if an array contains an index we must add 1 to it.
  32 C     EG: jjnr points to particles starting at 0
  33 C         type is indexed from 1 to ...
  34 C
  35
  36       subroutine FORLJC(ix,iy,iz,qi,
  37      $     pos,nj,type,jjnr,charge,nbfp,
  38      $     faction,fip,
  39      $     Vc,Vnb)
  40
  41       implicit none
  42
  43       real      ix,iy,iz,qi
  44       real      pos(*),charge(*),faction(*),fip(3)
  45       integer*4 nj,jjnr(*),type(*)
  46       real      Vc,Vnb,nbfp(*)
  47
  48       integer   k,jnr,j3,tj
  49       real      twelve,six
  50       real      fX,fY,fZ
  51       real      rijX,rijY,rijZ
  52       real      fijscal,vijcoul
  53       real      vctot,vnbtot
  54       real      rinv1,rinv2,rinv6
  55       real      fjx,fjy,fjz
  56       real      tx,ty,tz,vnb6,vnb12
  57
  58       parameter(twelve=12.0,six=6.0)
  59
  60       fX     = 0
  61       fY     = 0
  62       fZ     = 0
  63       vctot  = 0
  64       vnbtot = 0
  65
  66 cray compiler directive ignore vector dependencies
  67 c$dir ivdep
  68       do k=1,nj
  69          jnr   = jjnr(k)+1
  70          j3    = 3*jnr-2
  71          rijX  = ix - pos(j3)
  72          rijY  = iy - pos(j3+1)
  73          rijZ  = iz - pos(j3+2)
  74
  75          rinv1       = 1.0/sqrt((rijX*rijX)+(rijY*rijY)+(rijZ*rijZ))
  76          rinv2       = rinv1*rinv1
  77          rinv6       = rinv2*rinv2*rinv2
  78
  79          tj          = 2*type(jnr)+1
  80          vnb6        = nbfp(tj)*rinv6
  81          vnb12       = nbfp(tj+1)*rinv6*rinv6
  82          vijcoul     = qi*charge(jnr)*rinv1
  83
  84          vctot       = vctot+vijcoul
  85          vnbtot      = vnbtot+vnb12-vnb6
  86          fijscal     = (twelve*vnb12-six*vnb6+vijcoul)*rinv2
  87
  88          fjx           = faction(j3)
  89          tx            = rijX*fijscal
  90          fX            = fX + tx
  91          faction(j3)   = fjx - tx
  92          fjy           = faction(j3+1)
  93          ty            = rijY*fijscal
  94          fY            = fY + ty
  95          faction(j3+1) = fjy - ty
  96          fjz           = faction(j3+2)
  97          tz            = rijZ*fijscal
  98          fZ            = fZ + tz
  99          faction(j3+2) = fjz - tz
 100
 101       end do
 102
 103       fip(1) = fX
 104       fip(2) = fY
 105       fip(3) = fZ
 106       Vc     = Vc  + vctot
 107       Vnb    = Vnb + vnbtot
 108
 109       return
 110
 111       end
 112
 113       program main
 114
 115       implicit none
 116
 117       integer maxatom,maxx,maxlist,maxtype
 118       parameter(maxatom=1000,maxx=3*maxatom,maxlist=100)
 119       parameter(maxtype=19)
 120
 121       real*4 ix,iy,iz,qi,pos(maxx),faction(maxx),fip(3)
 122       real*4 charge(maxatom),nbfp(2*maxtype),Vc,Vnb
 123       integer type(maxatom),jjnr(maxlist)
 124       integer i,i3,j
 125
 126 c     setup benchmark
 127       do i=1,maxtype
 128          nbfp(2*i-1) = 1e-3
 129          nbfp(2*i)   = 1e-6
 130       end do
 131
 132       type(1)=1
 133       do i=2,maxatom
 134          type(i)=1+mod(type(i-1)+91,maxtype)
 135       end do
 136
 137       do i=1,maxatom
 138          i3=3*(i-1)
 139          pos(i3+1) = i
 140          pos(i3+2) = i
 141          pos(i3+3) = 1
 142
 143          charge(i) = mod(i,2)-0.5
 144       end do
 145
 146       jjnr(1) = 13
 147       do i=2,maxlist
 148          jjnr(i) = mod(jjnr(i-1)+13,maxatom)
 149       end do
 150
 151       ix = 0.0
 152       iy = 0.0
 153       iz = 0.0
 154       qi = 1.0
 155       Vc = 0.0
 156       Vnb= 0.0
 157
 158 c     run it
 159
 160
 161       do j=1,100
 162          do i=1,maxatom
 163
 164             call FORLJC(ix,iy,iz,qi,
 165      $           pos,maxlist,type,jjnr,charge,nbfp,
 166      $           faction,fip,
 167      $           Vc,Vnb)
 168
 169          end do
 170       end do
 171
 172       print *,Vc, Vnb
 173
 174       stop
 175       end