From 4846d8283e529542e2be1214a3e074d2e93eebef Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 25 Aug 2021 18:38:08 -0400 Subject: [PATCH] convert chain.f to fortran 90+ style free format file chain.f90 this is to maintain compatibility with some newer fortran compilers that do not support legacy style fortran by default anymore. --- cmake/Modules/Tools.cmake | 2 +- doc/src/Tools.rst | 2 +- tools/chain.f | 321 ------------------------------------ tools/chain.f90 | 333 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 335 insertions(+), 323 deletions(-) delete mode 100644 tools/chain.f create mode 100644 tools/chain.f90 diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index 29b21dc390..40d3048dcc 100644 --- a/cmake/Modules/Tools.cmake +++ b/cmake/Modules/Tools.cmake @@ -9,7 +9,7 @@ if(BUILD_TOOLS) check_language(Fortran) if(CMAKE_Fortran_COMPILER) enable_language(Fortran) - add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f) + add_executable(chain.x ${LAMMPS_TOOLS_DIR}/chain.f90) target_link_libraries(chain.x PRIVATE ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES}) install(TARGETS chain.x DESTINATION ${CMAKE_INSTALL_BINDIR}) else() diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index 72cfc7db9e..5822286e5f 100644 --- a/doc/src/Tools.rst +++ b/doc/src/Tools.rst @@ -172,7 +172,7 @@ Chris Lorenz (chris.lorenz at kcl.ac.uk), King's College London. chain tool ---------------------- -The file chain.f creates a LAMMPS data file containing bead-spring +The file chain.f90 creates a LAMMPS data file containing bead-spring polymer chains and/or monomer solvent atoms. It uses a text file containing chain definition parameters as an input. The created chains and solvent atoms can strongly overlap, so LAMMPS needs to run diff --git a/tools/chain.f b/tools/chain.f deleted file mode 100644 index 52481f5d0d..0000000000 --- a/tools/chain.f +++ /dev/null @@ -1,321 +0,0 @@ -c Create LAMMPS data file for collection of -c polymer bead-spring chains of various lengths and bead sizes -c Syntax: chain < def.chain > data.file -c def.chain is input file that specifies the chains -c data.file is output file that will be input for LAMMPS -c includes image flags in data file so chains can be unraveled later - - program chain - - integer swaptype - integer, allocatable :: nchain(:),nmonomer(:) - integer, allocatable :: ntype(:),nbondtype(:) - integer, allocatable :: type(:),molecule(:) - integer, allocatable :: imagex(:),imagey(:),imagez(:) - real*8, allocatable :: x(:),y(:),z(:) - real*8, allocatable :: bondlength(:),restrict(:) - common xprd,yprd,zprd,xboundlo,xboundhi, - $ yboundlo,yboundhi,zboundlo,zboundhi - real*8 random - 900 format(a) - 901 format(2f15.6,a) - 902 format(i3,f5.1) - 903 format(i10,i8,i8,3f10.4,3i4) - 904 format(i9,i3,2i9) - -c read chain definitions - - read (5,*) - read (5,*) - read (5,*) rhostar - read (5,*) iseed - read (5,*) nsets - read (5,*) swaptype - - allocate(nchain(nsets)) - allocate(nmonomer(nsets)) - allocate(ntype(nsets)) - allocate(nbondtype(nsets)) - allocate(bondlength(nsets)) - allocate(restrict(nsets)) - - do iset = 1,nsets - read (5,*) - read (5,*) nchain(iset) - read (5,*) nmonomer(iset) - read (5,*) ntype(iset) - read (5,*) nbondtype(iset) - read (5,*) bondlength(iset) - read (5,*) restrict(iset) - enddo - -c natoms = total # of monomers - - natoms = 0 - do iset = 1,nsets - natoms = natoms + nchain(iset)*nmonomer(iset) - enddo - - allocate(x(natoms)) - allocate(y(natoms)) - allocate(z(natoms)) - allocate(type(natoms)) - allocate(molecule(natoms)) - allocate(imagex(natoms)) - allocate(imagey(natoms)) - allocate(imagez(natoms)) - -c setup box size (sigma = 1.0) - - volume = natoms/rhostar - xprd = volume**(1.0/3.0) - yprd = xprd - zprd = xprd - - xboundlo = -xprd/2.0 - xboundhi = -xboundlo - yboundlo = xboundlo - yboundhi = xboundhi - zboundlo = xboundlo - zboundhi = xboundhi - -c generate random chains -c loop over sets and chains in each set - - n = 0 - nmolecule = 0 - - do iset = 1,nsets - do ichain = 1,nchain(iset) - nmolecule = nmolecule + 1 - -c random starting point for the chain in the box - - x1 = 0.0 - y1 = 0.0 - z1 = 0.0 - x2 = xboundlo + random(iseed)*xprd - y2 = yboundlo + random(iseed)*yprd - z2 = zboundlo + random(iseed)*zprd - -c store 1st monomer of chain -c 1st monomer is always in original box (image = 0) - - call pbc(x2,y2,z2) - n = n + 1 - x(n) = x2 - y(n) = y2 - z(n) = z2 - type(n) = ntype(iset) - imagex(n) = 0 - imagey(n) = 0 - imagez(n) = 0 - if (swaptype == 0) then - molecule(n) = nmolecule - else - molecule(n) = 1 - endif - -c generate rest of monomers in this chain - - do imonomer = 2,nmonomer(iset) - - x0 = x1 - y0 = y1 - z0 = z1 - x1 = x2 - y1 = y2 - z1 = z2 - -c random point inside sphere of unit radius - - 10 xinner = 2.0*random(iseed) - 1.0 - yinner = 2.0*random(iseed) - 1.0 - zinner = 2.0*random(iseed) - 1.0 - rsq = xinner*xinner + yinner*yinner + zinner*zinner - if (rsq > 1.0) goto 10 - -c project point to surface of sphere of unit radius - - r = sqrt(rsq) - xsurf = xinner/r - ysurf = yinner/r - zsurf = zinner/r - -c create new point by scaling unit offsets by bondlength (sigma = 1.0) - - x2 = x1 + xsurf*bondlength(iset) - y2 = y1 + ysurf*bondlength(iset) - z2 = z1 + zsurf*bondlength(iset) - -c check that new point meets restriction requirement -c only for 3rd monomer and beyond - - dx = x2 - x0 - dy = y2 - y0 - dz = z2 - z0 - r = sqrt(dx*dx + dy*dy + dz*dz) - - if (imonomer > 2 .and. r <= restrict(iset)) goto 10 - -c store new point -c if delta to previous bead is large, then increment/decrement image flag - - call pbc(x2,y2,z2) - n = n + 1 - x(n) = x2 - y(n) = y2 - z(n) = z2 - type(n) = ntype(iset) - - if (abs(x(n)-x(n-1)) < 2.0*bondlength(iset)) then - imagex(n) = imagex(n-1) - else if (x(n) - x(n-1) < 0.0) then - imagex(n) = imagex(n-1) + 1 - else if (x(n) - x(n-1) > 0.0) then - imagex(n) = imagex(n-1) - 1 - endif - - if (abs(y(n)-y(n-1)) < 2.0*bondlength(iset)) then - imagey(n) = imagey(n-1) - else if (y(n) - y(n-1) < 0.0) then - imagey(n) = imagey(n-1) + 1 - else if (y(n) - y(n-1) > 0.0) then - imagey(n) = imagey(n-1) - 1 - endif - - if (abs(z(n)-z(n-1)) < 2.0*bondlength(iset)) then - imagez(n) = imagez(n-1) - else if (z(n) - z(n-1) < 0.0) then - imagez(n) = imagez(n-1) + 1 - else if (z(n) - z(n-1) > 0.0) then - imagez(n) = imagez(n-1) - 1 - endif - - if (swaptype == 0) then - molecule(n) = nmolecule - else if (swaptype == 1) then - molecule(n) = imonomer - else if (swaptype == 2) then - if (imonomer <= nmonomer(iset)/2) then - molecule(n) = imonomer - else - molecule(n) = nmonomer(iset)+1-imonomer - endif - endif - - enddo - - enddo - enddo - -c compute quantities needed for LAMMPS file - - nbonds = 0 - ntypes = 0 - nbondtypes = 0 - do iset = 1,nsets - nbonds = nbonds + nchain(iset)*(nmonomer(iset)-1) - if (ntype(iset) > ntypes) ntypes = ntype(iset) - if (nbondtype(iset) > nbondtypes) - $ nbondtypes = nbondtype(iset) - enddo - -c write out LAMMPS file - - write (6,900) 'LAMMPS FENE chain 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,901) xboundlo,xboundhi,' xlo xhi' - write (6,901) yboundlo,yboundhi,' ylo yhi' - write (6,901) zboundlo,zboundhi,' zlo zhi' - - write (6,*) - write (6,900) 'Masses' - write (6,*) - - do i = 1,ntypes - write (6,902) i,1.0 - enddo - - write (6,*) - write (6,900) 'Atoms' - write (6,*) - - do i = 1,natoms - write (6,903) i,molecule(i),type(i),x(i),y(i),z(i), - $ imagex(i),imagey(i),imagez(i) - enddo - - if (nbonds > 0) then - - write (6,*) - write (6,900) 'Bonds' - write (6,*) - - n = 0 - m = 0 - do iset = 1,nsets - do ichain = 1,nchain(iset) - do imonomer = 1,nmonomer(iset) - n = n + 1 - if (imonomer /= nmonomer(iset)) then - m = m + 1 - write (6,904) m,nbondtype(iset),n,n+1 - endif - enddo - enddo - enddo - - endif - - end - -c ************ -c Subroutines -c ************ - -c periodic boundary conditions - map atom back into periodic box - - subroutine pbc(x,y,z) - common xprd,yprd,zprd,xboundlo,xboundhi, - $ yboundlo,yboundhi,zboundlo,zboundhi - - if (x < xboundlo) x = x + xprd - if (x >= xboundhi) x = x - xprd - if (y < yboundlo) y = y + yprd - if (y >= yboundhi) y = y - yprd - if (z < zboundlo) z = z + zprd - if (z >= zboundhi) z = z - zprd - - return - end - - -c RNG from Numerical Recipes - - real*8 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/chain.f90 b/tools/chain.f90 new file mode 100644 index 0000000000..1efb421cae --- /dev/null +++ b/tools/chain.f90 @@ -0,0 +1,333 @@ +! Create LAMMPS data file for collection of +! polymer bead-spring chains of various lengths and bead sizes +! Syntax: chain < def.chain > data.file +! def.chain is input file that specifies the chains +! data.file is output file that will be input for LAMMPS +! includes image flags in data file so chains can be unraveled later + +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,z) + REAL(KIND=8), INTENT(inout) :: x,y,z + + IF (x < xboundlo) x = x + xprd + IF (x >= xboundhi) x = x - xprd + IF (y < yboundlo) y = y + yprd + IF (y >= yboundhi) y = y - yprd + IF (z < zboundlo) z = z + zprd + IF (z >= zboundhi) z = z - zprd + + 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 chain + USE box + USE rng + IMPLICIT NONE + + INTEGER, ALLOCATABLE :: nchain(:),nmonomer(:) + INTEGER, ALLOCATABLE :: ntype(:),nbondtype(:) + INTEGER, ALLOCATABLE :: atomtype(:),molecule(:) + INTEGER, ALLOCATABLE :: imagex(:),imagey(:),imagez(:) + REAL(kind=8), ALLOCATABLE :: x(:),y(:),z(:) + REAL(kind=8), ALLOCATABLE :: bondlength(:),restrict(:) + INTEGER :: i, n, m, nmolecule, natoms, ntypes, nbonds, nbondtypes + INTEGER :: swaptype, iseed, nsets, iset, ichain, imonomer + REAL(kind=8) :: r, rhostar, volume, rsq, xinner, yinner, zinner, xsurf, ysurf, zsurf + REAL(kind=8) :: x0, y0, z0, x1, y1, z1, x2, y2, z2, dx, dy, dz + + LOGICAL :: again + + ! read chain definitions + + READ (5,*) + READ (5,*) + READ (5,*) rhostar + READ (5,*) iseed + READ (5,*) nsets + READ (5,*) swaptype + + ALLOCATE(nchain(nsets)) + ALLOCATE(nmonomer(nsets)) + ALLOCATE(ntype(nsets)) + ALLOCATE(nbondtype(nsets)) + ALLOCATE(bondlength(nsets)) + ALLOCATE(restrict(nsets)) + + DO iset = 1,nsets + READ (5,*) + READ (5,*) nchain(iset) + READ (5,*) nmonomer(iset) + READ (5,*) ntype(iset) + READ (5,*) nbondtype(iset) + READ (5,*) bondlength(iset) + READ (5,*) restrict(iset) + ENDDO + + ! natoms = total # of monomers + + natoms = 0 + DO iset = 1,nsets + natoms = natoms + nchain(iset)*nmonomer(iset) + ENDDO + + ALLOCATE(x(natoms)) + ALLOCATE(y(natoms)) + ALLOCATE(z(natoms)) + ALLOCATE(atomtype(natoms)) + ALLOCATE(molecule(natoms)) + ALLOCATE(imagex(natoms)) + ALLOCATE(imagey(natoms)) + ALLOCATE(imagez(natoms)) + + ! setup box size (sigma = 1.0) + + volume = natoms/rhostar + xprd = volume**(1.0/3.0) + yprd = xprd + zprd = xprd + + xboundlo = -xprd/2.0 + xboundhi = -xboundlo + yboundlo = xboundlo + yboundhi = xboundhi + zboundlo = xboundlo + zboundhi = xboundhi + + ! generate random chains + ! loop over sets and chains in each set + + n = 0 + nmolecule = 0 + + DO iset = 1,nsets + DO ichain = 1,nchain(iset) + nmolecule = nmolecule + 1 + + ! random starting point for the chain in the box + + x1 = 0.0 + y1 = 0.0 + z1 = 0.0 + x2 = xboundlo + random(iseed)*xprd + y2 = yboundlo + random(iseed)*yprd + z2 = zboundlo + random(iseed)*zprd + + ! store 1st monomer of chain + ! 1st monomer is always in original box (image = 0) + + CALL pbc(x2,y2,z2) + + n = n + 1 + x(n) = x2 + y(n) = y2 + z(n) = z2 + atomtype(n) = ntype(iset) + imagex(n) = 0 + imagey(n) = 0 + imagez(n) = 0 + IF (swaptype == 0) THEN + molecule(n) = nmolecule + ELSE + molecule(n) = 1 + END IF + + ! generate rest of monomers in this chain + + DO imonomer = 2, nmonomer(iset) + + x0 = x1 + y0 = y1 + z0 = z1 + x1 = x2 + y1 = y2 + z1 = z2 + + again = .TRUE. + DO WHILE (again) + ! random point inside sphere of unit radius + + xinner = 2.0*random(iseed) - 1.0 + yinner = 2.0*random(iseed) - 1.0 + zinner = 2.0*random(iseed) - 1.0 + rsq = xinner*xinner + yinner*yinner + zinner*zinner + IF (rsq > 1.0) CYCLE + + ! project point to surface of sphere of unit radius + + r = SQRT(rsq) + xsurf = xinner/r + ysurf = yinner/r + zsurf = zinner/r + + ! create new point by scaling unit offsets by bondlength (sigma = 1.0) + + x2 = x1 + xsurf*bondlength(iset) + y2 = y1 + ysurf*bondlength(iset) + z2 = z1 + zsurf*bondlength(iset) + + ! check that new point meets restriction requirement + ! only for 3rd monomer and beyond + + dx = x2 - x0 + dy = y2 - y0 + dz = z2 - z0 + r = SQRT(dx*dx + dy*dy + dz*dz) + + IF (imonomer > 2 .AND. r <= restrict(iset)) CYCLE + + ! store new point + again = .FALSE. + + ! if delta to previous bead is large, then increment/decrement image flag + + CALL pbc(x2,y2,z2) + n = n + 1 + x(n) = x2 + y(n) = y2 + z(n) = z2 + atomtype(n) = ntype(iset) + + IF (ABS(x(n)-x(n-1)) < 2.0*bondlength(iset)) THEN + imagex(n) = imagex(n-1) + ELSE IF (x(n) - x(n-1) < 0.0) THEN + imagex(n) = imagex(n-1) + 1 + ELSE IF (x(n) - x(n-1) > 0.0) THEN + imagex(n) = imagex(n-1) - 1 + ENDIF + + IF (ABS(y(n)-y(n-1)) < 2.0*bondlength(iset)) THEN + imagey(n) = imagey(n-1) + ELSE IF (y(n) - y(n-1) < 0.0) THEN + imagey(n) = imagey(n-1) + 1 + ELSE IF (y(n) - y(n-1) > 0.0) THEN + imagey(n) = imagey(n-1) - 1 + ENDIF + + IF (ABS(z(n)-z(n-1)) < 2.0*bondlength(iset)) THEN + imagez(n) = imagez(n-1) + ELSE IF (z(n) - z(n-1) < 0.0) THEN + imagez(n) = imagez(n-1) + 1 + ELSE IF (z(n) - z(n-1) > 0.0) THEN + imagez(n) = imagez(n-1) - 1 + ENDIF + + IF (swaptype == 0) THEN + molecule(n) = nmolecule + ELSE IF (swaptype == 1) THEN + molecule(n) = imonomer + ELSE IF (swaptype == 2) THEN + IF (imonomer <= nmonomer(iset)/2) THEN + molecule(n) = imonomer + ELSE + molecule(n) = nmonomer(iset)+1-imonomer + ENDIF + ENDIF + ENDDO + ENDDO + + ENDDO + ENDDO + + ! compute quantities needed for LAMMPS file + + nbonds = 0 + ntypes = 0 + nbondtypes = 0 + DO iset = 1,nsets + nbonds = nbonds + nchain(iset)*(nmonomer(iset)-1) + IF (ntype(iset) > ntypes) ntypes = ntype(iset) + IF (nbondtype(iset) > nbondtypes) nbondtypes = nbondtype(iset) + ENDDO + + ! write out LAMMPS file + + WRITE (6,*) 'LAMMPS FENE chain 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,'(2F15.6,A)') xboundlo,xboundhi,' xlo xhi' + WRITE (6,'(2F15.6,A)') yboundlo,yboundhi,' ylo yhi' + WRITE (6,'(2F15.6,A)') zboundlo,zboundhi,' zlo zhi' + + WRITE (6,*) + WRITE (6,*) 'Masses' + WRITE (6,*) + + DO i = 1,ntypes + WRITE (6,'(i3,f5.1)') i,1.0 + ENDDO + + WRITE (6,*) + WRITE (6,*) 'Atoms' + WRITE (6,*) + + DO i = 1,natoms + WRITE (6,'(I10,I8,I8,3F10.4,3I4)') i,molecule(i),atomtype(i),x(i),y(i),z(i), & + imagex(i),imagey(i),imagez(i) + ENDDO + + IF (nbonds > 0) THEN + WRITE (6,*) + WRITE (6,*) 'Bonds' + WRITE (6,*) + + n = 0 + m = 0 + DO iset = 1,nsets + DO ichain = 1,nchain(iset) + DO imonomer = 1,nmonomer(iset) + n = n + 1 + IF (imonomer /= nmonomer(iset)) THEN + m = m + 1 + WRITE (6,'(i9,i3,2i9)') m,nbondtype(iset),n,n+1 + ENDIF + ENDDO + ENDDO + ENDDO + ENDIF + + DEALLOCATE(nchain, nmonomer, ntype, nbondtype, bondlength, restrict) + DEALLOCATE(x, y, z, atomtype, molecule, imagex, imagey, imagez) + +END PROGRAM chain