1
subroutine newcfg(natom1,natom2,itask,rcluster,rcut)
3
integer natom1, natom2, itask
4
integer n, nc, nx, ny, nz, ntot
6
parameter ( nc = 50, n = 2*nc ** 3)
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
20
c *******************************************************
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
42
100 format('md.cfg',i1)
43
101 format('md.cfg',i2)
44
102 format('md.cfg',i3)
45
103 format('md.cfg',i4)
47
open (unit=2,file=filename,status='unknown',form='formatted')
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
65
multx = mult * float(nside)
66
multy = mult * float(nside)
67
multz = mult * float(nside)
70
dx = multx/float(nside)
71
dy = multy/float(nside)
72
dz = multz/float(nside)
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
85
iy = int(sqrt(dble(irad)**2-dble(ix-ixcent)**2))
89
iz = int(sqrt(dble(irad)**2-dble(ix-ixcent)**2
90
+ -dble(iy-iycent)**2))
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
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)
104
if (fill(ix,iy,iz,2).eq.0.and.icnt.lt.natom2) then
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)
114
if (icnt.lt.natom2) then
119
c Find center of mass of cluster
125
rx = rx + coords(natom1+icnt,1)
126
ry = ry + coords(natom1+icnt,2)
127
rz = rz + coords(natom1+icnt,3)
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)
134
c Find radius of cluster
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
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
155
coords(icnt,1) = xadd + lj(1,1)
156
coords(icnt,2) = yadd + lj(2,1)
157
coords(icnt,3) = zadd + lj(3,1)
160
if (fill(ix,iy,iz,2).eq.0.and.icnt.lt.natom1) then
162
coords(icnt,1) = xadd + lj(1,2)
163
coords(icnt,2) = yadd + lj(2,2)
164
coords(icnt,3) = zadd + lj(3,2)
173
c If there are no atoms of type 1, then assume cluster is
174
c isolated and set size to 100 LJ units
176
if (natom1.eq.0) then
181
if (rcluster.lt.rcut) then
184
write(2,2100) natom1+natom2
185
write(2,2000) multx,multy,multz,rcluster
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)
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)
199
200 format (i5,3f12.5)