      PROGRAM genetic
      IMPLICIT NONE
      REAL*8 Eselect(1000)
      INTEGER ntot, L, L2, T(400,2), i, igeneration, iswap,
     +       ngeneration, nmutation, npop, nrecombine
      CHARACTER*80 s(1000)
      CALL initial(s,L,L2,T,npop,nrecombine,nmutation,ngeneration)
      DO igeneration = 1, ngeneration
         ntot = npop
         DO iswap = 1, nrecombine
            CALL recombine(s,L2,npop,ntot)
         END DO
         DO i = 1, nmutation
            CALL mutate(s,L2,npop,ntot)
         END DO
         CALL selection(s,L,L2,T,npop,ntot,Eselect)
      END DO
      CALL showoutput(s,L2,npop,Eselect)
      END
C
      SUBROUTINE initial(s,L,L2,T,npop,nrecombine,nmutation
     +,ngeneration)
      IMPLICIT NONE
      REAL dmy, rnd
      INTEGER L, L2, T(400,2), npop, nrecombine, nmutation,
     +        ngeneration, i, ipop, j
      CHARACTER*80 s(1000)
      dmy = rnd(-1)
      WRITE(*,'(A,$)') 'string (lattice) size = '
      READ(*,*) L
      L2 = L*L
      WRITE(*,'(A,$)') 'population number = '
      READ(*,*) npop
*     L is linear dimension of phenotype
      WRITE(*,'(A,$)') 'number of recombinations per generation = '
      READ(*,*) nrecombine
      WRITE(*,'(A,$)') 'number of mutations per generation = '
      READ(*,*) nmutation
      WRITE(*,'(A,$)') 'number of generations = '
      READ(*,*) ngeneration
*     create random population of genotypes
      DO ipop = 1, npop
         s(ipop) = ''
         DO i = 1, L2
            IF(rnd(0) .GT. 0.5D0) THEN
               s(ipop)(i:i) = '1'
            ELSE
               s(ipop)(i:i) = '0'
            END IF
         END DO
      END DO
*     create random bond network of Tij's
      DO i = 1, L2
         DO j = 1, 2
            IF(rnd(0) .GT. 0.5D0) THEN
               T(i,j) = 1
            ELSE
               T(i,j) = -1
            END IF
         END DO
      END DO
      END
C
      SUBROUTINE recombine(s,L2,npop,ntot)
*     choose two strings (genotypes) to recombine
      IMPLICIT NONE
      REAL rnd
      INTEGER pbc, L2, npop, i, j, pos, size, ntot
      CHARACTER*80 s1, s2, s(1000)
      LOGICAL Ltmp1_
      i = int(npop*rnd(0)) + 1
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         j = int(npop*rnd(0)) + 1
         Ltmp1_ = i .EQ. j
      END DO
      size = int(rnd(0)*L2/2) + 1
      pos = int(rnd(0)*L2) + 1
      s1 = s(i)
      s2 = s(j)
      IF(pos + size .LE. L2) THEN
         s(ntot+1) = s1(1:pos-1)//s2(pos:pos+size)//s1(pos+size+1:L2)
         s(ntot+2) = s2(1:pos-1)//s1(pos:pos+size)//s2(pos+size+1:L2)
      ELSE
         pbc = pos + size - L2
         s(ntot+1) = s2(1:pbc)//s1(pbc+1:pos-1)//s2(pos:L2)
         s(ntot+2) = s1(1:pbc)//s2(pbc+1:pos-1)//s1(pos:L2)
      END IF
      ntot = ntot + 2
      END
C
      SUBROUTINE mutate(s,L2,npop,ntot)
      IMPLICIT NONE
      REAL rnd
      INTEGER L2, npop, i, pos, ntot
      CHARACTER*80 c*1, s(1000)
      i = int(rnd(0)*npop) + 1
      pos = int(rnd(0)*L2) + 1
      c = s(i)(pos:pos)
      IF(c .EQ. '1') THEN
         c = '0'
      ELSE
         c = '1'
      END IF
      s(ntot + 1) = s(i)(1:pos-1)//c//s(i)(pos+1:L2)
      ntot = ntot + 1
      END
C
      SUBROUTINE convert(a,L,ss)
*     coverts strings of 0's and 1's to 2D array of spins
*     that is, genotype to phenotype
      IMPLICIT NONE
      INTEGER n, L, i, j, ss(30,30)
      CHARACTER*80 a
      DO i = 1, L
         DO j = 1, L
            n = (j-1)*L + i
            IF(a(n:n) .EQ. '1') THEN
               ss(i,j) = 1
            ELSE
               ss(i,j) = -1
            END IF
         END DO
      END DO
      END
C
      SUBROUTINE energy(L,ss,T,E)
      IMPLICIT NONE
      REAL*8 E
      INTEGER L, T(400,2), i, ip, j, jp, n,  ss(30,30)
      E = 0
      DO i = 1, L
         ip = i + 1
         IF(ip .GT. L) ip = 1
         DO j = 1, L
            jp = j + 1
            IF(jp .GT. L) jp = 1
            n = (j-1)*L + i
            E = E - T(n,1)*ss(i,j)*ss(ip,j) - T(n,2)*ss(i,j)*ss(i,jp)
         END DO
      END DO
      END
C
      SUBROUTINE selection(s,L,L2,T,npop,ntot,Eselect)
      IMPLICIT NONE
      REAL*8 E, Elist(0:1000), Esum, Eselect(1000)
      REAL rnd
      INTEGER L, L2, T(400,2), npop, i, ipop, ntot, ss(30,30)
      CHARACTER*80 s(1000), sv(1000)
      LOGICAL Ltmp1_
      Esum = 0
      Elist(0) = 0
      DO i = 1, ntot
         CALL convert(s(i),L,ss)
         CALL energy(L,ss,T,E)
         Esum = Esum -  E + 2*L2            ! contribution to Esum > 0
         Elist(i) = Esum
      END DO
      DO i = 1, ntot
         sv(i) = s(i)
      END DO
*     select new population
      DO ipop = 1, npop
         E = Esum*rnd(0)
         i = 0
         Ltmp1_ = .TRUE.
         DO WHILE(Ltmp1_)
            i = i + 1
            Ltmp1_ = E .GE. Elist(i)
         END DO
         s(ipop) = sv(i)
         Eselect(ipop) = Elist(i-1) - Elist(i) + 2*L2
      END DO
      END
C
      SUBROUTINE showoutput(s,L2,npop,Eselect)
      IMPLICIT NONE
      REAL*8 Eselect(1000)
      INTEGER npop, ipop, L2
      CHARACTER*80 s(1000)
      DO ipop = 1, npop
         WRITE(*,'(I6,3X,A,F12.6)') ipop, s(ipop)(1:L2), Eselect(ipop)
      END DO
      WRITE(*,*)
      END
