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