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