From 0aded3931b16d79288c692851ea6938b9a58d3b8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 26 Aug 2021 00:36:03 -0400 Subject: [PATCH] convert to Fortran 2003 --- tools/micelle2d.f | 233 -------------------------------------------- tools/micelle2d.f90 | 231 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 231 insertions(+), 233 deletions(-) delete mode 100644 tools/micelle2d.f create mode 100644 tools/micelle2d.f90 diff --git a/tools/micelle2d.f b/tools/micelle2d.f deleted file mode 100644 index 7bc879fcf9..0000000000 --- a/tools/micelle2d.f +++ /dev/null @@ -1,233 +0,0 @@ -c Create LAMMPS 2003 input file for 2d LJ simulation of micelles -c -c Syntax: micelle2d < def.micelle2d > data.file -c -c def file contains size of system and number to turn into surfactants -c attaches a random surfactant tail(s) to random head -c if nonflag is set, will attach 2nd-neighbor bonds in surfactant -c solvent atoms = type 1 -c micelle heads = type 2 -c micelle tails = type 3,4,5,etc. - - program micelle2d - parameter (maxatom=100000,maxbond=10000) - real*4 x(2,maxatom) - integer type(maxatom),molecule(maxatom) - integer bondatom(2,maxbond),bondtype(maxbond) - common xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi, - $ zboundlo,zboundhi - 999 format (3i7,3f8.3) - 998 format (4i7) - - read (5,*) - read (5,*) - read (5,*) rhostar - read (5,*) iseed - read (5,*) nx,ny - read (5,*) nsurf - read (5,*) r0 - read (5,*) ntails - read (5,*) nonflag - - natoms = nx*ny - if (natoms+nsurf*ntails.gt.maxatom) then - write (6,*) 'Too many atoms - boost maxatom' - call exit(0) - endif - - nbonds = nsurf*ntails - if (nonflag.eq.1) nbonds = nbonds + nsurf*(ntails-1) - if (nbonds.gt.maxbond) then - write (6,*) 'Too many surfactants - boost maxbond' - call exit(0) - endif - -c box size - - rlattice = (1.0/rhostar) ** 0.5 - - xboundlo = 0.0 - xboundhi = nx*rlattice - yboundlo = 0.0 - yboundhi = ny*rlattice - zboundlo = -0.1 - zboundhi = 0.1 - - sigma = 1.0 - - xprd = xboundhi - xboundlo - yprd = yboundhi - yboundlo - zprd = zboundhi - zboundlo - -c initial square lattice of solvents - - m = 0 - do j = 1,ny - do i = 1,nx - m = m + 1 - x(1,m) = xboundlo + (i-1)*rlattice - x(2,m) = yboundlo + (j-1)*rlattice - molecule(m) = 0 - type(m) = 1 - enddo - enddo - -c turn some into surfactants with molecule ID -c head changes to type 2 -c create ntails for each head of types 3,4,... -c each tail is at distance r0 away in straight line with random orientation - - do i = 1,nsurf - - 10 m = random(iseed)*natoms + 1 - if (m.gt.natoms) m = natoms - if (molecule(m) .ne. 0) goto 10 - - molecule(m) = i - type(m) = 2 - - angle = random(iseed)*2.0*3.1415926 - do j = 1,ntails - k = (i-1)*ntails + j - x(1,natoms+k) = x(1,m) + cos(angle)*j*r0*sigma - x(2,natoms+k) = x(2,m) + sin(angle)*j*r0*sigma - molecule(natoms+k) = i - type(natoms+k) = 2+j - call pbc(x(1,natoms+k),x(2,natoms+k)) - if (j.eq.1) bondatom(1,k) = m - if (j.ne.1) bondatom(1,k) = natoms+k-1 - bondatom(2,k) = natoms+k - bondtype(k) = 1 - enddo - - enddo - -c if nonflag is set, add (ntails-1) 2nd nearest neighbor bonds to end -c of bond list -c k = location in bondatom list where nearest neighbor bonds for -c this surfactant are stored - - if (nonflag.eq.1) then - - nbonds = nsurf*ntails - do i = 1,nsurf - do j = 1,ntails-1 - k = (i-1)*ntails + j - nbonds = nbonds + 1 - bondatom(1,nbonds) = bondatom(1,k) - bondatom(2,nbonds) = bondatom(2,k+1) - bondtype(nbonds) = 2 - enddo - enddo - - endif - -c write LAMMPS data file - - natoms = natoms + nsurf*ntails - nbonds = nsurf*ntails - if (nonflag.eq.1) nbonds = nbonds + nsurf*(ntails-1) - ntypes = 2 + ntails - nbondtypes = 1 - if (nonflag.eq.1) nbondtypes = 2 - - if (nsurf.eq.0) then - ntypes = 1 - nbondtypes = 0 - endif - - write (6,*) 'LAMMPS 2d micelle data file' - write (6,*) - - write (6,*) natoms,' atoms' - write (6,*) nbonds,' bonds' - write (6,*) 0,' angles' - write (6,*) 0,' dihedrals' - write (6,*) 0,' impropers' - write (6,*) - - write (6,*) ntypes,' atom types' - write (6,*) nbondtypes,' bond types' - write (6,*) 0,' angle types' - write (6,*) 0,' dihedral types' - write (6,*) 0,' improper types' - write (6,*) - - write (6,*) xboundlo,xboundhi,' xlo xhi' - write (6,*) yboundlo,yboundhi,' ylo yhi' - write (6,*) zboundlo,zboundhi,' zlo zhi' - - write (6,*) - write (6,*) 'Masses' - write (6,*) - - do i = 1,ntypes - write (6,*) i,1.0 - enddo - - write (6,*) - write (6,*) 'Atoms' - write (6,*) - - do i = 1,natoms - write (6,999) i,molecule(i),type(i),x(1,i),x(2,i),0.0 - enddo - - if (nsurf.gt.0) then - - write (6,*) - write (6,*) 'Bonds' - write (6,*) - - do i = 1,nbonds - write (6,998) i,bondtype(i),bondatom(1,i),bondatom(2,i) - enddo - - endif - -c write Xmovie bond geometry file - - open(1,file='bond.micelle') - - write (1,*) 'ITEM: BONDS' - do i = 1,nbonds - write (1,*) bondtype(i),bondatom(1,i),bondatom(2,i) - enddo - - close(1) - - end - - -c ************ -c Subroutines -c ************ - -c periodic boundary conditions - map atom back into periodic box - - subroutine pbc(x,y) - common xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi, - $ zboundlo,zboundhi - - if (x.lt.xboundlo) x = x + xprd - if (x.ge.xboundhi) x = x - xprd - if (y.lt.yboundlo) y = y + yprd - if (y.ge.yboundhi) y = y - yprd - - return - end - - -c RNG - compute in double precision, return single - - real*4 function random(iseed) - real*8 aa,mm,sseed - parameter (aa=16807.0D0,mm=2147483647.0D0) - - sseed = iseed - sseed = mod(aa*sseed,mm) - random = sseed/mm - iseed = sseed - - return - end diff --git a/tools/micelle2d.f90 b/tools/micelle2d.f90 new file mode 100644 index 0000000000..19d8e513c7 --- /dev/null +++ b/tools/micelle2d.f90 @@ -0,0 +1,231 @@ +! Create LAMMPS data file for 2d LJ simulation of micelles +! +! Syntax: micelle2d < def.micelle2d > data.file +! +! def file contains size of system and number to turn into surfactants +! attaches a random surfactant tail(s) to random head +! if nonflag is set, will attach 2nd-neighbor bonds in surfactant +! solvent atoms = type 1 +! micelle heads = type 2 +! micelle tails = type 3,4,5,etc. + +MODULE box + IMPLICIT NONE + PUBLIC + REAL(KIND=8) :: xprd,yprd,zprd,xboundlo,xboundhi,yboundlo,yboundhi,zboundlo,zboundhi + +CONTAINS + + ! periodic boundary conditions - map atom back into periodic box + + SUBROUTINE pbc(x,y) + REAL(KIND=8), INTENT(inout) :: x,y + + IF (x < xboundlo) x = x + xprd + IF (x >= xboundhi) x = x - xprd + IF (y < yboundlo) y = y + yprd + IF (y >= yboundhi) y = y - yprd + + END SUBROUTINE pbc +END MODULE box + +MODULE rng + IMPLICIT NONE + +CONTAINS + + ! *very* minimal random number generator + + REAL(KIND=8) FUNCTION random(iseed) + IMPLICIT NONE + INTEGER, INTENT(inout) :: iseed + REAL(KIND=8), PARAMETER :: aa=16807.0_8, mm=2147483647.0_8 + REAL(KIND=8) :: sseed + + sseed = REAL(iseed) + sseed = MOD(aa*sseed,mm) + random = sseed/mm + iseed = INT(sseed) + END FUNCTION random +END MODULE rng + +PROGRAM micelle2d + USE box + USE rng + IMPLICIT NONE + + REAL(kind=8), ALLOCATABLE :: x(:,:) + INTEGER, ALLOCATABLE :: atomtype(:), molecule(:) + INTEGER, ALLOCATABLE :: bondatom(:,:),bondtype(:) + INTEGER :: natoms, maxatom, ntypes, nbonds, nbondtypes, iseed + INTEGER :: i, j, k, m, nx, ny, nsurf, ntails, nonflag + REAL(kind=8) :: rhostar, rlattice, sigma, angle,r0 + REAL(kind=8), parameter :: pi = 3.14159265358979323846_8 + LOGICAL :: again + + READ(5,*) + READ(5,*) + READ(5,*) rhostar + READ(5,*) iseed + READ(5,*) nx,ny + READ(5,*) nsurf + READ(5,*) r0 + READ(5,*) ntails + READ(5,*) nonflag + + natoms = nx*ny + maxatom = natoms + nsurf*ntails + ALLOCATE(x(2,maxatom), molecule(maxatom), atomtype(maxatom)) + + nbonds = nsurf*ntails + IF (nonflag.EQ.1) nbonds = nbonds + nsurf*(ntails-1) + ALLOCATE(bondatom(2,nbonds), bondtype(nbonds)) + +! box size + + rlattice = (1.0_8/rhostar) ** 0.5_8 + + xboundlo = 0.0_8 + xboundhi = nx*rlattice + yboundlo = 0.0_8 + yboundhi = ny*rlattice + zboundlo = -0.1_8 + zboundhi = 0.1_8 + + sigma = 1.0_8 + + xprd = xboundhi - xboundlo + yprd = yboundhi - yboundlo + zprd = zboundhi - zboundlo + +! initial square lattice of solvents + + m = 0 + DO j = 1,ny + DO i = 1,nx + m = m + 1 + x(1,m) = xboundlo + (i-1)*rlattice + x(2,m) = yboundlo + (j-1)*rlattice + molecule(m) = 0 + atomtype(m) = 1 + ENDDO + ENDDO + +! turn some into surfactants with molecule ID +! head changes to type 2 +! create ntails for each head of types 3,4,... +! each tail is at distance r0 away in straight line with random orientation + + DO i = 1,nsurf + + again = .TRUE. + DO WHILE(again) + m = INT(random(iseed)*natoms + 1) + IF (m > natoms) m = natoms + IF (molecule(m) /= 0) CYCLE + again = .FALSE. + END DO + molecule(m) = i + atomtype(m) = 2 + + angle = random(iseed)*2.0_8*pi + DO j = 1,ntails + k = (i-1)*ntails + j + x(1,natoms+k) = x(1,m) + COS(angle)*j*r0*sigma + x(2,natoms+k) = x(2,m) + SIN(angle)*j*r0*sigma + molecule(natoms+k) = i + atomtype(natoms+k) = 2+j + CALL pbc(x(1,natoms+k),x(2,natoms+k)) + IF (j == 1) bondatom(1,k) = m + IF (j /= 1) bondatom(1,k) = natoms+k-1 + bondatom(2,k) = natoms+k + bondtype(k) = 1 + ENDDO + + ENDDO + +! if nonflag is set, add (ntails-1) 2nd nearest neighbor bonds to end +! of bond list +! k = location in bondatom list where nearest neighbor bonds for +! this surfactant are stored + + IF (nonflag == 1) THEN + + nbonds = nsurf*ntails + DO i = 1,nsurf + DO j = 1,ntails-1 + k = (i-1)*ntails + j + nbonds = nbonds + 1 + bondatom(1,nbonds) = bondatom(1,k) + bondatom(2,nbonds) = bondatom(2,k+1) + bondtype(nbonds) = 2 + ENDDO + ENDDO + + ENDIF + +! write LAMMPS data file + + natoms = natoms + nsurf*ntails + nbonds = nsurf*ntails + IF (nonflag == 1) nbonds = nbonds + nsurf*(ntails-1) + ntypes = 2 + ntails + nbondtypes = 1 + IF (nonflag == 1) nbondtypes = 2 + + IF (nsurf == 0) THEN + ntypes = 1 + nbondtypes = 0 + ENDIF + + WRITE (6,*) 'LAMMPS 2d micelle data file' + WRITE (6,*) + + WRITE (6,*) natoms,' atoms' + WRITE (6,*) nbonds,' bonds' + WRITE (6,*) 0,' angles' + WRITE (6,*) 0,' dihedrals' + WRITE (6,*) 0,' impropers' + WRITE (6,*) + + WRITE (6,*) ntypes,' atom types' + WRITE (6,*) nbondtypes,' bond types' + WRITE (6,*) 0,' angle types' + WRITE (6,*) 0,' dihedral types' + WRITE (6,*) 0,' improper types' + WRITE (6,*) + + WRITE (6,*) xboundlo,xboundhi,' xlo xhi' + WRITE (6,*) yboundlo,yboundhi,' ylo yhi' + WRITE (6,*) zboundlo,zboundhi,' zlo zhi' + + WRITE (6,*) + WRITE (6,*) 'Masses' + WRITE (6,*) + + DO i = 1,ntypes + WRITE (6,*) i,1.0 + ENDDO + + WRITE (6,*) + WRITE (6,*) 'Atoms # molecular' + WRITE (6,*) + + DO i = 1,natoms + WRITE (6,'(3I7,3F8.3)') i,molecule(i),atomtype(i),x(1,i),x(2,i),0.0 + ENDDO + + IF (nsurf > 0) THEN + + WRITE (6,*) + WRITE (6,*) 'Bonds' + WRITE (6,*) + + DO i = 1,nbonds + WRITE (6,'(4I7)') i,bondtype(i),bondatom(1,i),bondatom(2,i) + ENDDO + + ENDIF + + DEALLOCATE(x,molecule,atomtype,bondtype,bondatom) +END PROGRAM micelle2d