1! Create LAMMPS data file for collection of 2! polymer bead-spring chains of various lengths and bead sizes 3! Syntax: chain < def.chain > data.file 4! def.chain is input file that specifies the chains 5! data.file is output file that will be input for LAMMPS 6! includes image flags in data file so chains can be unraveled later 7 8MODULE boxchain 9 IMPLICIT NONE 10 PUBLIC 11 REAL(KIND=8) :: xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,zboundlo,zboundhi 12 13CONTAINS 14 15 ! periodic boundary conditions - map atom back into periodic box 16 17 SUBROUTINE pbc(x,y,z) 18 REAL(KIND=8), INTENT(inout) :: x,y,z 19 20 IF (x < xboundlo) x = x + xprd 21 IF (x >= xboundhi) x = x - xprd 22 IF (y < yboundlo) y = y + yprd 23 IF (y >= yboundhi) y = y - yprd 24 IF (z < zboundlo) z = z + zprd 25 IF (z >= zboundhi) z = z - zprd 26 27 END SUBROUTINE pbc 28END MODULE boxchain 29 30MODULE rngchain 31 IMPLICIT NONE 32 33CONTAINS 34 35 ! *very* minimal random number generator 36 37 REAL(KIND=8) FUNCTION random(iseed) 38 IMPLICIT NONE 39 INTEGER, INTENT(inout) :: iseed 40 REAL(KIND=8), PARAMETER :: aa=16807.0_8, mm=2147483647.0_8 41 REAL(KIND=8) :: sseed 42 43 sseed = REAL(iseed) 44 sseed = MOD(aa*sseed,mm) 45 random = sseed/mm 46 iseed = INT(sseed) 47 END FUNCTION random 48END MODULE rngchain 49 50PROGRAM chain 51 USE boxchain 52 USE rngchain 53 IMPLICIT NONE 54 55 INTEGER, ALLOCATABLE :: nchain(:),nmonomer(:) 56 INTEGER, ALLOCATABLE :: ntype(:),nbondtype(:) 57 INTEGER, ALLOCATABLE :: atomtype(:),molecule(:) 58 INTEGER, ALLOCATABLE :: imagex(:),imagey(:),imagez(:) 59 REAL(KIND=8), ALLOCATABLE :: x(:),y(:),z(:) 60 REAL(KIND=8), ALLOCATABLE :: bondlength(:),restrict(:) 61 INTEGER :: i, n, m, nmolecule, natoms, ntypes, nbonds, nbondtypes 62 INTEGER :: swaptype, iseed, nsets, iset, ichain, imonomer 63 REAL(KIND=8) :: r, rhostar, volume, rsq, xinner, yinner, zinner, xsurf, ysurf, zsurf 64 REAL(KIND=8) :: x0, y0, z0, x1, y1, z1, x2, y2, z2, dx, dy, dz 65 66 LOGICAL :: again 67 68 ! read chain definitions 69 70 READ (5,*) 71 READ (5,*) 72 READ (5,*) rhostar 73 READ (5,*) iseed 74 READ (5,*) nsets 75 READ (5,*) swaptype 76 77 ALLOCATE(nchain(nsets)) 78 ALLOCATE(nmonomer(nsets)) 79 ALLOCATE(ntype(nsets)) 80 ALLOCATE(nbondtype(nsets)) 81 ALLOCATE(bondlength(nsets)) 82 ALLOCATE(restrict(nsets)) 83 84 DO iset = 1,nsets 85 READ (5,*) 86 READ (5,*) nchain(iset) 87 READ (5,*) nmonomer(iset) 88 READ (5,*) ntype(iset) 89 READ (5,*) nbondtype(iset) 90 READ (5,*) bondlength(iset) 91 READ (5,*) restrict(iset) 92 ENDDO 93 94 ! natoms = total # of monomers 95 96 natoms = 0 97 DO iset = 1,nsets 98 natoms = natoms + nchain(iset)*nmonomer(iset) 99 ENDDO 100 101 ALLOCATE(x(natoms)) 102 ALLOCATE(y(natoms)) 103 ALLOCATE(z(natoms)) 104 ALLOCATE(atomtype(natoms)) 105 ALLOCATE(molecule(natoms)) 106 ALLOCATE(imagex(natoms)) 107 ALLOCATE(imagey(natoms)) 108 ALLOCATE(imagez(natoms)) 109 110 ! setup box size (sigma = 1.0) 111 112 volume = natoms/rhostar 113 xprd = volume**(1.0/3.0) 114 yprd = xprd 115 zprd = xprd 116 117 xboundlo = -xprd/2.0 118 xboundhi = -xboundlo 119 yboundlo = xboundlo 120 yboundhi = xboundhi 121 zboundlo = xboundlo 122 zboundhi = xboundhi 123 124 ! generate random chains 125 ! loop over sets and chains in each set 126 127 n = 0 128 nmolecule = 0 129 130 DO iset = 1,nsets 131 DO ichain = 1,nchain(iset) 132 nmolecule = nmolecule + 1 133 134 ! random starting point for the chain in the box 135 136 x1 = 0.0 137 y1 = 0.0 138 z1 = 0.0 139 x2 = xboundlo + random(iseed)*xprd 140 y2 = yboundlo + random(iseed)*yprd 141 z2 = zboundlo + random(iseed)*zprd 142 143 ! store 1st monomer of chain 144 ! 1st monomer is always in original box (image = 0) 145 146 CALL pbc(x2,y2,z2) 147 148 n = n + 1 149 x(n) = x2 150 y(n) = y2 151 z(n) = z2 152 atomtype(n) = ntype(iset) 153 imagex(n) = 0 154 imagey(n) = 0 155 imagez(n) = 0 156 IF (swaptype == 0) THEN 157 molecule(n) = nmolecule 158 ELSE 159 molecule(n) = 1 160 END IF 161 162 ! generate rest of monomers in this chain 163 164 DO imonomer = 2, nmonomer(iset) 165 166 x0 = x1 167 y0 = y1 168 z0 = z1 169 x1 = x2 170 y1 = y2 171 z1 = z2 172 173 again = .TRUE. 174 DO WHILE (again) 175 ! random point inside sphere of unit radius 176 177 xinner = 2.0*random(iseed) - 1.0 178 yinner = 2.0*random(iseed) - 1.0 179 zinner = 2.0*random(iseed) - 1.0 180 rsq = xinner*xinner + yinner*yinner + zinner*zinner 181 IF (rsq > 1.0) CYCLE 182 183 ! project point to surface of sphere of unit radius 184 185 r = SQRT(rsq) 186 xsurf = xinner/r 187 ysurf = yinner/r 188 zsurf = zinner/r 189 190 ! create new point by scaling unit offsets by bondlength (sigma = 1.0) 191 192 x2 = x1 + xsurf*bondlength(iset) 193 y2 = y1 + ysurf*bondlength(iset) 194 z2 = z1 + zsurf*bondlength(iset) 195 196 ! check that new point meets restriction requirement 197 ! only for 3rd monomer and beyond 198 199 dx = x2 - x0 200 dy = y2 - y0 201 dz = z2 - z0 202 r = SQRT(dx*dx + dy*dy + dz*dz) 203 204 IF (imonomer > 2 .AND. r <= restrict(iset)) CYCLE 205 206 ! store new point 207 again = .FALSE. 208 209 ! if delta to previous bead is large, then increment/decrement image flag 210 211 CALL pbc(x2,y2,z2) 212 n = n + 1 213 x(n) = x2 214 y(n) = y2 215 z(n) = z2 216 atomtype(n) = ntype(iset) 217 218 IF (ABS(x(n)-x(n-1)) < 2.0*bondlength(iset)) THEN 219 imagex(n) = imagex(n-1) 220 ELSE IF (x(n) - x(n-1) < 0.0) THEN 221 imagex(n) = imagex(n-1) + 1 222 ELSE IF (x(n) - x(n-1) > 0.0) THEN 223 imagex(n) = imagex(n-1) - 1 224 ENDIF 225 226 IF (ABS(y(n)-y(n-1)) < 2.0*bondlength(iset)) THEN 227 imagey(n) = imagey(n-1) 228 ELSE IF (y(n) - y(n-1) < 0.0) THEN 229 imagey(n) = imagey(n-1) + 1 230 ELSE IF (y(n) - y(n-1) > 0.0) THEN 231 imagey(n) = imagey(n-1) - 1 232 ENDIF 233 234 IF (ABS(z(n)-z(n-1)) < 2.0*bondlength(iset)) THEN 235 imagez(n) = imagez(n-1) 236 ELSE IF (z(n) - z(n-1) < 0.0) THEN 237 imagez(n) = imagez(n-1) + 1 238 ELSE IF (z(n) - z(n-1) > 0.0) THEN 239 imagez(n) = imagez(n-1) - 1 240 ENDIF 241 242 IF (swaptype == 0) THEN 243 molecule(n) = nmolecule 244 ELSE IF (swaptype == 1) THEN 245 molecule(n) = imonomer 246 ELSE IF (swaptype == 2) THEN 247 IF (imonomer <= nmonomer(iset)/2) THEN 248 molecule(n) = imonomer 249 ELSE 250 molecule(n) = nmonomer(iset)+1-imonomer 251 ENDIF 252 ENDIF 253 ENDDO 254 ENDDO 255 256 ENDDO 257 ENDDO 258 259 ! compute quantities needed for LAMMPS file 260 261 nbonds = 0 262 ntypes = 0 263 nbondtypes = 0 264 DO iset = 1,nsets 265 nbonds = nbonds + nchain(iset)*(nmonomer(iset)-1) 266 IF (ntype(iset) > ntypes) ntypes = ntype(iset) 267 IF (nbondtype(iset) > nbondtypes) nbondtypes = nbondtype(iset) 268 ENDDO 269 270 ! write out LAMMPS file 271 272 WRITE (6,*) 'LAMMPS FENE chain data file' 273 WRITE (6,*) 274 275 WRITE (6,*) natoms,' atoms' 276 WRITE (6,*) nbonds,' bonds' 277 WRITE (6,*) 0,' angles' 278 WRITE (6,*) 0,' dihedrals' 279 WRITE (6,*) 0,' impropers' 280 WRITE (6,*) 281 282 WRITE (6,*) ntypes,' atom types' 283 WRITE (6,*) nbondtypes,' bond types' 284 WRITE (6,*) 0,' angle types' 285 WRITE (6,*) 0,' dihedral types' 286 WRITE (6,*) 0,' improper types' 287 WRITE (6,*) 288 289 WRITE (6,'(2F15.6,A)') xboundlo,xboundhi,' xlo xhi' 290 WRITE (6,'(2F15.6,A)') yboundlo,yboundhi,' ylo yhi' 291 WRITE (6,'(2F15.6,A)') zboundlo,zboundhi,' zlo zhi' 292 293 WRITE (6,*) 294 WRITE (6,*) 'Masses' 295 WRITE (6,*) 296 297 DO i = 1,ntypes 298 WRITE (6,'(i3,f5.1)') i,1.0 299 ENDDO 300 301 WRITE (6,*) 302 WRITE (6,*) 'Atoms # molecular' 303 WRITE (6,*) 304 305 DO i = 1,natoms 306 WRITE (6,'(I10,I8,I8,3F10.4,3I4)') i,molecule(i),atomtype(i),x(i),y(i),z(i), & 307 imagex(i),imagey(i),imagez(i) 308 ENDDO 309 310 IF (nbonds > 0) THEN 311 WRITE (6,*) 312 WRITE (6,*) 'Bonds' 313 WRITE (6,*) 314 315 n = 0 316 m = 0 317 DO iset = 1,nsets 318 DO ichain = 1,nchain(iset) 319 DO imonomer = 1,nmonomer(iset) 320 n = n + 1 321 IF (imonomer /= nmonomer(iset)) THEN 322 m = m + 1 323 WRITE (6,'(i9,i3,2i9)') m,nbondtype(iset),n,n+1 324 ENDIF 325 ENDDO 326 ENDDO 327 ENDDO 328 ENDIF 329 330 DEALLOCATE(nchain, nmonomer, ntype, nbondtype, bondlength, restrict) 331 DEALLOCATE(x, y, z, atomtype, molecule, imagex, imagey, imagez) 332 333END PROGRAM chain 334