~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/tools/ga-4-3/global/examples/md_cluster/newcfg.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      subroutine newcfg(natom1,natom2,itask,rcluster,rcut)
 
2
      implicit none 
 
3
      integer natom1, natom2, itask
 
4
      integer     n, nc, nx, ny, nz, ntot
 
5
 
 
6
      parameter ( nc = 50, n = 2*nc ** 3)
 
7
 
 
8
      double precision        lj(3,2)
 
9
      double precision        cell,cell2,mult,xadd,yadd,zadd
 
10
      double precision        multx,multy,multz,cmx,cmy,cmz
 
11
      double precision        x,y,z,dx,dy,dz,rcluster,r,rx,ry,rz,rmax
 
12
      double precision        volume, radius, pi, rcut
 
13
      double precision        coords(n, 3)
 
14
      integer     i, j, ix, iy, iz, iref, m, mp, mp1, mp4
 
15
      integer     one,two,nside, irad, ixcent, iycent, izcent
 
16
      integer     ixmin,ixmax,iymin,iymax,izmin,izmax
 
17
      integer     fill(nc,nc,nc,2), icnt
 
18
      character*32 filename
 
19
c
 
20
c   *******************************************************
 
21
c
 
22
      mult = 1.5d00
 
23
      lj(1,1) = 0.0d00
 
24
      lj(2,1) = 0.0d00
 
25
      lj(3,1) = 0.0d00
 
26
      lj(1,2) = mult/2.0d00
 
27
      lj(2,2) = mult/2.0d00
 
28
      lj(3,2) = mult/2.0d00
 
29
c
 
30
      one = 1
 
31
      two = 2
 
32
c
 
33
      if (itask.lt.10) then
 
34
        write(filename,100) itask
 
35
      else if (itask.ge.10.and.itask.lt.100) then
 
36
        write(filename,101) itask
 
37
      else if (itask.ge.100.and.itask.lt.1000) then
 
38
        write(filename,102) itask
 
39
      else if (itask.ge.1000.and.itask.lt.10000) then
 
40
        write(filename,103) itask
 
41
      endif
 
42
  100 format('md.cfg',i1)
 
43
  101 format('md.cfg',i2)
 
44
  102 format('md.cfg',i3)
 
45
  103 format('md.cfg',i4)
 
46
 
 
47
      open (unit=2,file=filename,status='unknown',form='formatted')
 
48
c
 
49
      nside = nint((0.5d00*dble(natom1+natom2))**(1.0d00/3.0d00))
 
50
      if (2*nside**3.lt.natom1+natom2) nside = nside+1
 
51
      pi = 4.0d00*atan(1.0d00)
 
52
      volume = 0.5d00*dble(natom2)*mult**3
 
53
      radius = 3.0d00*volume/(4.0d00*pi)
 
54
      radius = radius**(1.0d00/3.0d00)
 
55
      irad = int(radius/mult) + 1
 
56
      if (2*irad.ge.nside) nside = 2*irad+2
 
57
      do ix = 1, nside
 
58
        do iy = 1, nside
 
59
          do iz = 1, nside
 
60
            fill(ix,iy,iz,1) = 0
 
61
            fill(ix,iy,iz,2) = 0
 
62
          end do
 
63
        end do
 
64
      end do
 
65
      multx = mult * float(nside)
 
66
      multy = mult * float(nside)
 
67
      multz = mult * float(nside)
 
68
      cell = mult
 
69
      cell2 = 0.5d00 * cell
 
70
      dx = multx/float(nside)
 
71
      dy = multy/float(nside)
 
72
      dz = multz/float(nside)
 
73
 
 
74
      volume = 0.5d00*dble(natom2)*mult**3
 
75
      radius = 3.0d00*volume/(4.0d00*pi)
 
76
      radius = radius**(1.0d00/3.0d00)
 
77
      irad = int(radius/mult) + 1
 
78
      ixcent = nside/2
 
79
      iycent = nside/2
 
80
      izcent = nside/2
 
81
      ixmin = ixcent - irad
 
82
      ixmax = ixcent + irad
 
83
 300  icnt = 0
 
84
      do ix = ixmin, ixmax
 
85
        iy = int(sqrt(dble(irad)**2-dble(ix-ixcent)**2))
 
86
        iymin = iycent - iy
 
87
        iymax = iycent + iy
 
88
        do iy = iymin, iymax
 
89
          iz = int(sqrt(dble(irad)**2-dble(ix-ixcent)**2
 
90
     +                               -dble(iy-iycent)**2))
 
91
          izmin = izcent - iz
 
92
          izmax = izcent + iz
 
93
          do iz = izmin, izmax
 
94
            xadd = dx*float(ix-1) - multx/2.0d00
 
95
            yadd = dy*float(iy-1) - multy/2.0d00
 
96
            zadd = dz*float(iz-1) - multz/2.0d00
 
97
            if (fill(ix,iy,iz,1).eq.0.and.icnt.lt.natom2) then
 
98
              icnt = icnt + 1
 
99
              coords(natom1+icnt,1) = xadd + lj(1,1)
 
100
              coords(natom1+icnt,2) = yadd + lj(2,1)
 
101
              coords(natom1+icnt,3) = zadd + lj(3,1)
 
102
              fill(ix,iy,iz,1) = 1
 
103
            endif
 
104
            if (fill(ix,iy,iz,2).eq.0.and.icnt.lt.natom2) then
 
105
              icnt = icnt + 1
 
106
              coords(natom1+icnt,1) = xadd + lj(1,2)
 
107
              coords(natom1+icnt,2) = yadd + lj(2,2)
 
108
              coords(natom1+icnt,3) = zadd + lj(3,2)
 
109
              fill(ix,iy,iz,2) = 1
 
110
            endif
 
111
          end do
 
112
        end do
 
113
      end do
 
114
      if (icnt.lt.natom2) then
 
115
        irad = irad+1
 
116
        go to 300
 
117
      endif
 
118
c
 
119
c  Find center of mass of cluster
 
120
c
 
121
      rx = 0.0d00
 
122
      ry = 0.0d00
 
123
      rz = 0.0d00
 
124
      do icnt = 1, natom2
 
125
        rx = rx + coords(natom1+icnt,1)
 
126
        ry = ry + coords(natom1+icnt,2)
 
127
        rz = rz + coords(natom1+icnt,3)
 
128
      end do
 
129
      cmx = rx/dble(natom2)
 
130
      cmy = ry/dble(natom2)
 
131
      cmz = rz/dble(natom2)
 
132
      rmax = rmax + sqrt(dx**2+dy**2+dz**2)
 
133
c
 
134
c  Find radius of cluster
 
135
c
 
136
      rmax = 0.0d00
 
137
      do icnt = 1, natom2
 
138
        rx = coords(natom1+icnt,1) - cmx
 
139
        ry = coords(natom1+icnt,2) - cmy
 
140
        rz = coords(natom1+icnt,3) - cmz
 
141
        r = sqrt(rx**2 + ry**2 + rz**2)
 
142
        if (r.gt.rmax) rmax = r
 
143
      end do
 
144
      rmax = rmax + 0.1
 
145
c
 
146
      icnt = 0
 
147
      do ix = 1, nside
 
148
        do iy = 1, nside
 
149
          do iz = 1, nside
 
150
            xadd = dx*float(ix-1) - multx/2.0d00
 
151
            yadd = dy*float(iy-1) - multy/2.0d00
 
152
            zadd = dz*float(iz-1) - multz/2.0d00
 
153
            if (fill(ix,iy,iz,1).eq.0.and.icnt.lt.natom1) then
 
154
              icnt = icnt + 1
 
155
              coords(icnt,1) = xadd + lj(1,1)
 
156
              coords(icnt,2) = yadd + lj(2,1)
 
157
              coords(icnt,3) = zadd + lj(3,1)
 
158
              fill(ix,iy,iz,1) = 1
 
159
            endif
 
160
            if (fill(ix,iy,iz,2).eq.0.and.icnt.lt.natom1) then
 
161
              icnt = icnt + 1
 
162
              coords(icnt,1) = xadd + lj(1,2)
 
163
              coords(icnt,2) = yadd + lj(2,2)
 
164
              coords(icnt,3) = zadd + lj(3,2)
 
165
              fill(ix,iy,iz,2) = 1
 
166
            endif
 
167
          end do
 
168
        end do
 
169
      end do
 
170
c
 
171
      rcluster = rmax
 
172
c
 
173
c  If there are no atoms of type 1, then assume cluster is
 
174
c  isolated and set size to 100 LJ units
 
175
c
 
176
      if (natom1.eq.0) then
 
177
        multx = 25.0
 
178
        multy = 25.0
 
179
        multz = 25.0
 
180
      endif
 
181
      if (rcluster.lt.rcut) then
 
182
        rcluster = rcut
 
183
      endif
 
184
      write(2,2100) natom1+natom2
 
185
      write(2,2000) multx,multy,multz,rcluster
 
186
      do i = 1, natom1
 
187
        coords(i,1) = coords(i,1) - cmx
 
188
        coords(i,2) = coords(i,2) - cmy
 
189
        coords(i,3) = coords(i,3) - cmz
 
190
        write(2,200) one,(coords(i,j),j=1,3)
 
191
      end do
 
192
      do i = natom1+1, natom1+natom2
 
193
        coords(i,1) = coords(i,1) - cmx
 
194
        coords(i,2) = coords(i,2) - cmy
 
195
        coords(i,3) = coords(i,3) - cmz
 
196
        write(2,200) two,(coords(i,j),j=1,3)
 
197
      end do
 
198
c
 
199
 200  format (i5,3f12.5)
 
200
2000  format (4f16.8)
 
201
2100  format (i8)
 
202
      close(2)
 
203
 
 
204
      return
 
205
      end