ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

convert chain.f to fortran 90+ style free format file chain.f90

Browse Source 
this is to maintain compatibility with some newer fortran compilers
that do not support legacy style fortran by default anymore.
This commit is contained in:
Axel Kohlmeyer
2021-08-25 18:38:08 -04:00
parent abd4a6cfa3
commit 4846d8283e
4 changed files with 335 additions and 323 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/Tools.cmake
View File

@ -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()

2
doc/src/Tools.rst
View File

@ -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

321
tools/chain.f
View File

@ -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

333
tools/chain.f90 Normal file
View File

@ -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