program gamatch

c     Program to index TPMCl triclinic crystal faces by using 
c     standard GAs to maximise the 1/deviation of the measured & 
c     observed interplanar angle 

c     Declaration
c
c     popsize - popoulation size
c     pc - probability of crossover
c     pm - probability of mutation at current generation
c     pmi - initial mutation probability
c     dpm - increase in mutation probability throughout the evolution
c     m - length of chromosome
c     index - position where the chromosome split to different faces
c     v - chromosome for the current generation
c     fv - fitness value for the current generation
c     vp - chromosome for the next generation
c     ncount - generation number
c     mcount - max. generation number
c     error - accepted fitness value
c     bigf - total population fitness value
c     rwf, rwp - probability vector/matrix
c     nv - no. of faces to represent the chromosome
c     pos - position to crossover (2 pts.)
c     pe - probability to inject the best old string to the population
c     prf - no. of faces to be specified
c     pr - position (start from the beginning of the chromosome) of the 
c          string to be protected from any change (= prf * 5)
c     face - MIs assigened to the faces
c     prnt - print flag (=1, print immediate results in crt)
c     disk - write flag (=1, save the immediate results in disk)
c     rmsflag - use rmsd% (spoil) = 1 or rmsd 

      integer popsize, m, mcount, nv
      parameter(popsize=200, m=15, mcount=3000, nv=3)

      integer v(popsize,m),vp(popsize,m),mv(m),mg,i,j,k,l,
     _        ncount,cp(popsize),s(m),ind(20),nd,pr,disk,
     _        plane(50,2),lookup(32,4),x(50,3),pos(2),lendir,
     _        face(10,3),face1(3),binno(5),prf,prnt,rmsflag

      double precision pc,pm,fv(popsize),mfv,error,pmi,dpm,
     _                 rwf(popsize),rwp(m),ps(popsize),z,pe,
     _                 qi(popsize),dummy,obsang(50),phi(50),
     _                 as,bs,cs,alps,bets,gams,scale,penality,
     _                 spoil,rmsd

      character*40 dire
      character*55 fnamein, fnameout

c     Assign operational parameters

      pc = 0.1
      pmi = 0.15
      pm = pmi
      dpm = 0.05
      pe = 0.1
      scale = 100.0
      penality = 0.8
      error = scale/0.4
      rmsflag = 0

c     Assign file names

      dire = '/home/joule/rgc/tam/GAs/'
      lendir = 24
      fnamein(1:lendir) = dire
      fnameout(1:lendir) = dire

      fnamein(lendir+1:lendir+11) = '2xtbfdh.dat'
      fnameout(lendir+1:lendir+6) = '2a.dat'

c     Assign print/write flag

      prnt = 1
      disk = 1

c     Assign protected face(s)

      prf = 0

      face(1,1) = 0 
      face(1,2) = 1
      face(1,3) = 0

      face(2,1) = -1
      face(2,2) = 1
      face(2,3) = 0

      face(3,1) = 1
      face(3,2) = 0
      face(3,3) = 0

      pr = 5 * prf

c     Get data file of observed angles & faces

      open(unit=13, file=fnamein, status='old')

      j = 1
10    read(13, *, err=20) plane(j,1),plane(j,2),obsang(j)
          j = j + 1
          goto 10
20    close(unit = 13)
      nd = j - 1

c     Create data file to record the evolution process

      open(unit=14, file=fnameout, status='unknown')

c     Set up the lookup table for difference faces

      do i = 0, 31
         lookup(i+1, 1) = i
      end do
      lookup(1, 2) = -1 
      lookup(1, 3) = 1
      lookup(1, 4) = 1
      lookup(2, 2) = 1 
      lookup(2, 3) = -1
      lookup(2, 4) = -1
      lookup(3, 2) = -1 
      lookup(3, 3) = 1
      lookup(3, 4) = 0
      lookup(4, 2) = 1 
      lookup(4, 3) = -1
      lookup(4, 4) = 0
      lookup(5, 2) = 0 
      lookup(5, 3) = 1
      lookup(5, 4) = 0
      lookup(6, 2) = 0 
      lookup(6, 3) = -1
      lookup(6, 4) = 0
      lookup(7, 2) = 1 
      lookup(7, 3) = 0
      lookup(7, 4) = 0
      lookup(8, 2) = -1 
      lookup(8, 3) = 0
      lookup(8, 4) = 0
      lookup(9, 2) = 0 
      lookup(9, 3) = 1
      lookup(9, 4) = 0
      lookup(10, 2) = 0 
      lookup(10, 3) = -1
      lookup(10, 4) = 0
      lookup(11, 2) = 0 
      lookup(11, 3) = 0
      lookup(11, 4) = 1
      lookup(12, 2) = 0 
      lookup(12, 3) = 0
      lookup(12, 4) = -1
      lookup(13, 2) = 0 
      lookup(13, 3) = 1
      lookup(13, 4) = 1
      lookup(14, 2) = 0 
      lookup(14, 3) = -1
      lookup(14, 4) = -1
      lookup(15, 2) = 1 
      lookup(15, 3) = 0
      lookup(15, 4) = 1
      lookup(16, 2) = -1 
      lookup(16, 3) = 0
      lookup(16, 4) = -1
      lookup(17, 2) = 1 
      lookup(17, 3) = 1
      lookup(17, 4) = 0
      lookup(18, 2) = -1 
      lookup(18, 3) = -1
      lookup(18, 4) = 0
      lookup(19, 2) = 0 
      lookup(19, 3) = -1
      lookup(19, 4) = 1
      lookup(20, 2) = 0 
      lookup(20, 3) = 1
      lookup(20, 4) = -1
      lookup(21, 2) = -1 
      lookup(21, 3) = 0
      lookup(21, 4) = 1
      lookup(22, 2) = 1 
      lookup(22, 3) = 0
      lookup(22, 4) = -1
      lookup(23, 2) = -1 
      lookup(23, 3) = 1
      lookup(23, 4) = 0
      lookup(24, 2) = 1 
      lookup(24, 3) = -1
      lookup(24, 4) = 0
      lookup(25, 2) = 1 
      lookup(25, 3) = 1
      lookup(25, 4) = 1
      lookup(26, 2) = -1 
      lookup(26, 3) = -1
      lookup(26, 4) = -1
      lookup(27, 2) = -1 
      lookup(27, 3) = 1
      lookup(27, 4) = 1
      lookup(28, 2) = 1 
      lookup(28, 3) = -1
      lookup(28, 4) = -1
      lookup(29, 2) = -1 
      lookup(29, 3) = -1
      lookup(29, 4) = 1
      lookup(30, 2) = 1 
      lookup(30, 3) = 1
      lookup(30, 4) = -1
      lookup(31, 2) = -1 
      lookup(31, 3) = 1
      lookup(31, 4) = -1
      lookup(32, 2) = 1 
      lookup(32, 3) = -1
      lookup(32, 4) = 1

      j = 0
      do i = 5, 5*(nv - 1), 5
         j = j + 1
         ind(j) = i
      end do
      call calpar(as,bs,cs,alps,bets,gams)

c     Protect the gene(s)

      do j = 1, prf
         do k = 1, 3
            face1(k) = face(j, k)
         end do 
         call face2bin(lookup, face1,
     _                               binno)
         do i = 1, popsize
            do l = (j - 1) * 5 + 1, j * 5
               v(i, l) = binno(l - (j - 1) * 5)
               vp(i, l) = v(i, l)
            end do
         end do
      end do

c     Initialise the initial state of the random number generator

      call g05ccf

c     Establish initial population

      ncount = 0
      mg = 0
      do i = 1, popsize
         do j = pr + 1, m
             call random(v(i,j), dummy)
         end do
      end do  

c     Start evolution

      if (prnt.eq.1) then
         print*, '   Gen.  fitness  Total fitness'
      end if
      if (disk.eq.1) then
         write(14,*)'   Gen.  fitness  Total fitness'  
      end if

100   continue
          ncount = ncount + 1
          bigf = 0.0      
          do i = 1, popsize
              do j = 1, m
                 s(j) = v(i, j)
              end do
              call angle(s,m,ind,nv,nd,as,bs,cs,alps,bets,
     _              gams,obsang,plane,lookup,scale,penality,
     _              spoil,rmsd,x,phi)
              if (rmsflag.eq.1) then
                  fv(i) = spoil
                else
                  fv(i) = rmsd
              end if
              bigf = bigf + fv(i)
              call random(j, rwf(i))
          end do
          call fmax(popsize, fv, 
     _                          dummy, i)
          if ((ncount.eq.1).or.(dummy.gt.mfv)) then
             mfv = dummy
             do j = 1, m
                mv(j) = v(i,j)
             end do
             mg = ncount
          end if

c     Stochastic remainder selection

c       qi - scaled fitness
c       ps - remainder

          dummy = bigf/popsize
          j = 0
          do i = 1, popsize
             qi(i) = fv(i)/dummy
             pos(1) = int(qi(i))
             ps(i) = qi(i) - pos(1)
             do k = 1, pos(1)
                j = j + 1
                if (j.gt.popsize) then
                   goto 130
                end if
                do l = 1, m                 
                   vp(j, l) = v(i, l)
                end do
             end do
          end do
          do i = j, popsize
125          continue
             call random(k, dummy)
             cp(1) = int((popsize - 1) * dummy) + 1
             call random(k, dummy)
             if (dummy.lt.ps(cp(1))) then
                 do l = 1, m
                    vp(i, l) = v(cp(1), l)
                 end do
                 ps(cp(1)) = 0.0
               else
                 goto 125
             end if
          end do
                    
130       continue

c     Crossover

          j = 0
          do i = 1, popsize
             call random(k, rwf(i))
             if (rwf(i).lt.pc) then
                 j = j + 1
                 cp(j) = i
             end if 
          end do
          if (real(j/2).ne.(j/2.0)) then
             j = j + 1
             call random(k, dummy)
             cp(j) = int((popsize - 1) * dummy) + 1
          end if                
          do i = 2, j, 2
             do l = 1, 2
                call random(k, dummy)
                pos(l) = int((m - 1 - pr) * dummy) + 1 + pr
             end do

c            2 points crossover

             if (pos(1).lt.pos(2)) then
                do k = pos(1), pos(2)
                   s(k) = vp(cp(i - 1), k)
                   vp(cp(i - 1), k) = vp(cp(i), k)
                   vp(cp(i), k) = s(k)
                end do
              else
                do k = pr + 1, pos(2)
                   s(k) = vp(cp(i - 1), k)
                   vp(cp(i - 1), k) = vp(cp(i), k)
                   vp(cp(i), k) = s(k)
                end do
                do k = pos(1), m
                   s(k) = vp(cp(i - 1), k)
                   vp(cp(i - 1), k) = vp(cp(i), k)
                   vp(cp(i), k) = s(k)
                end do
             end if
          end do  

c     Elitism

c     Inject the best old chromosome into the population

          call random(k, dummy)
          if (dummy.lt.pe) then
            do i = 1, popsize
              do j = 1, m
                 s(j) = vp(i, j)
              end do
              call angle(s,m,ind,nv,nd,as,bs,cs,alps,bets,
     _            gams,obsang,plane,lookup,scale,penality,
     _            spoil,rmsd,x,phi)
              if (rmsflag.eq.1) then
                  fv(i) = spoil
               else
                  fv(i) = rmsd
              end if
            end do
            call fmax(popsize, fv, 
     _                          z, l)
            if (mfv.ge.z) then
               call fmin(popzise, fv, dummy, j)
               do k = 1, m
                  vp(j, k) = mv(k)
               end do
             else
               mfv = z
               do k = 1, m
                  mv(k) = vp(l, k)
               end do
               mg = ncount
            end if
          end if

c      Mutation

          pm = pm + dpm / mcount
          do i = 1, popsize
             l = 0
             do j = pr + 1, m             
                call random(k, rwp(j))
                if (rwp(j).lt.pm) then
                   l = l + 1
                   s(l) = j
                end if
             end do
             do j = 1, l
                call random(k, dummy)
                vp(i,s(j)) = k
             end do
          end do

c      Update the population 

          do i = 1, popsize
             do j = 1, m
                v(i, j) = vp(i, j)
             end do
          end do
          if (prnt.eq.1) then
              print 1000, ncount, mfv, bigf
          end if
          if (disk.eq.1) then
              write(14,1000) ncount, mfv, bigf 
          end if
          if ((ncount.lt.mcount).and.(mfv.lt.error)) then
              goto 100
          end if  
200   continue

      do i = 1, m
         s(i) = mv(i)
      end do
      call angle(s,m,ind,nv,nd,as,bs,cs,alps,bets,gams,
     _            obsang,plane,lookup,scale,penality,
     _            spoil,rmsd,x,phi)
      if (rmsflag.eq.1) then
          dummy = spoil
        else
          dummy = rmsd
      end if
      if (prnt.eq.1) then
         print*, ' '
         print*, 'Minimised angular deviation/fitness value'
         print 1010, scale/rmsd, dummy
         print*, 'ncount (optm)'
         print*, ' '
         print*, mg
         print*, '     Assigned faces'
      end if

      write(14,*) ' '
      write(14,*) 'Minimised angular deviation/fitness value'
      write(14,1010) scale/rmsd, dummy
      write(14,*) 'ncount (optm)'
      write(14,*) mg
      write(14,*) ' '
      write(14,*) '     Assigned faces'
      do i = 1, nv
            if (prnt.eq.1) then
               print 1020, i, x(i,1), x(i,2), x(i,3)
            end if
            write(14,1020) i, x(i,1), x(i,2), x(i,3)
      end do
      if (prnt.eq.1) then
         print*, ' '
         print*, '     Faces     Obs.    Calc.'
      end if
      write(14,*) ' '
      write(14,*) '     Faces     Obs.    Calc.'
      do i = 1, nd
            if (prnt.eq.1) then
               print 1030,plane(i,1),plane(i,2),obsang(i),phi(i) 
            end if
            write(14,1030) plane(i,1),plane(i,2),obsang(i),phi(i)
      end do 

      write(14,*) ' '
      write(14,1040) 'Population size', popsize 
      write(14,1040) 'Max. generation', mcount
      write(14,1050) 'Pr(crossover)  ', pc
      write(14,1050) 'Pr(mutation,0) ', pmi
      write(14,1050) 'd Pr(mutation) ', dpm
      write(14,1050) 'Pr(elitism)    ', pe
      write(14,1040) 'Protected face ', prf
      write(14,1050) 'Penalty        ', penality
      write(14,1060) 'Data file name ', fnamein
      write(14,1060) 'O/P file name  ', fnameout

      close(unit = 14)

1000  format(i5,2x,f10.2,2x,f10.2)
1010  format(f7.2,20x,f7.2)
1020  format(i5,i5,i5,i5)
1030  format(i5,i5,2x,f7.1,2x,f7.2)
1040  format(a20,2x,i7)
1050  format(a20,2x,f7.2)
1060  format(a20,a55)

      end

      subroutine angle(s,m,ind,nv,nd,as,bs,cs,alps,bets,gams,
     _                  obsang,plane,lookup,scale,penality,
     _                  spoil,rmsd,x,phi)

      integer s(1000),ind(50),m,ni,i,j,lookup(32,4),plane(50,2),
     _        ix(50),x(50,3),nv,nd
      integer bin2dec 

      double precision as,bs,cs,alps,bets,gams,spoil,h(3),hp(3),
     _            rmsd,phi(50),obsang(50),penality,factor,scale
      double precision iang,rmill

c     s - chromosome in binary format 
c     index - indicate the position to split s into different variables
c        e.g. [1,1,1,1,1,1,1,1,1,1] to 2 variables at position 5, ind = 5
c        [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1] at 5 & 10 then ind = [5, 10]

      ni = nv - 1
      factor = 1.0

      ix(1) = bin2dec(s, 1, ind(1))
      do i = 2, ni
         j = 1 + ind(i - 1)
         ix(i) = bin2dec(s, j, ind(i))
      end do
      j = 1 + ind(ni)
      ix(ni+1) = bin2dec(s, j, m)
      do i = 1, nv - 1
         do j = i + 1, nv
            if (ix(i).eq.ix(j)) then
                factor = factor * penality 
             else
                if ((ix(i).le.5).or.(ix(j).le.5)) then
                    call chkeqgene(ix(i),ix(j),0,26,penality,
     _                             factor)
                    call chkeqgene(ix(i),ix(j),1,27,penality,
     _                             factor)
                    call chkeqgene(ix(i),ix(j),2,22,penality,
     _                             factor)
                    call chkeqgene(ix(i),ix(j),3,23,penality,
     _                             factor)
                    call chkeqgene(ix(i),ix(j),4,8,penality,
     _                             factor)
                    call chkeqgene(ix(i),ix(j),5,9,penality,
     _                             factor)
                end if    
            end if
         end do
      end do
      do i = 1, nv
         do j = 1, 32
            if (ix(i).eq.lookup(j,1)) then
               x(i,1) = lookup(j,2)
               x(i,2) = lookup(j,3)
               x(i,3) = lookup(j,4)
            end if
         end do 
      end do

      spoil = 0.0
      rmsd = 0.0
      do i = 1, nd
         do j = 1, 3
           h(j)=rmill(x(plane(i,1),j))
           hp(j)=rmill(x(plane(i,2),j))
         end do
         phi(i) = iang(h,hp,as,bs,cs,alps,bets,gams)
         spoil = spoil + ((obsang(i) - phi(i))/obsang(i))**2
         rmsd = rmsd + (obsang(i) - phi(i))**2
      end do
      spoil = factor*scale/(sqrt(spoil/nd))
      rmsd = factor*scale/(sqrt(rmsd/nd))

      return
      end

      subroutine chkeqgene(ix1, ix2, a, b, penality,
     _                                             factor)
      
c     check the repeated genes: 0,1,0; -1,1,1; -1,1,0 & their complement

      integer ix1, ix2, a, b
      double precision penality, factor

      if (((ix1.eq.a).and.(ix2.eq.b)).or.((ix2.eq.a).
     _                                and.(ix1.eq.b))) then
         factor = factor * penality
      end if

      return
      end

      integer function bin2dec(s, i1, i2)

c     convert binary number to integer
c
c     s - i/p binary number
c     i1 - start bit position
c     i2 - end bit position

      integer i1, i2, i, dummy, s(1000)

      dummy=0
      do i = i1, i2
          dummy = dummy + s(i)*2**(i2-i)
      end do
      bin2dec = dummy

      return
      end

      subroutine face2bin(lookup,face,
     _                                binno)

c     convert the face to the corresponding binary value in the lookup table

      integer i, lookup(32, 4), face(3), binno(5), j

      do i = 1, 32
         if (((lookup(i,2).eq.face(1)).and.(lookup(i,3).eq.face(2))).
     _        and.(lookup(i,4).eq.face(3))) then
              j = lookup(i, 1)
              goto 2000
         end if
      end do
2000  if (i.eq.33) then
        do i = 1, 5
           binno(i) = -1
        end do
       else 
        do i = 0, 3   
           dummy = j / 2
           j = j - 2 * dummy    
           binno(5-i) = j 
           j = dummy
        end do 
        binno(1) = j
      end if 

      return
      end 

      subroutine random(numi, numd)

c     generate uniformly dirtributed random no. within 0 & 1
c     numi - integer o/p (0 or 1)
c     numd - double precision o/p

      integer numi
      double precision numd

      numd = 10.0*(g05daf(0.9,1.0)-0.9)
      if (numd.ge.0.5) then
           numi = 1
        else
           numi = 0
      end if

      return
      end

      subroutine fmax(np, y,
     _                       maxi, pos)

c     find max. in an array

      integer i, np, pos
      double precision y(10000), maxi
      
      maxi = y(1)
      pos = 1
      do i = 2, np
         if (y(i).gt.maxi) then
             maxi = y(i)
             pos = i
         end if
      end do

      return
      end

      subroutine fmin(np, y,
     _                       mini, i)

c     find min. in an array

      integer i, np, pos
      double precision y(10000), mini

      mini = y(1)
      pos = 1
      do i = 2, np
         if (y(i).lt.mini) then 
            mini = y(i)
            pos = i
         end if
      end do

      return
      end

      subroutine calpar(as,bs,cs,alps,bets,gams)

c     calculate reciprocal lattice parameters for triclinic TPMCl
c
c     formular taken from:
c
c     International Table For X-ray Crystallography, vol II,
c        Kynoch Press, Brmingham, 1973, p106.
c
c     Declaration of variables
c
c     a, b, c, alp, bet, gam - unit cell parameters
c     v - unit cell volume
c     s - (alp + bet + gam)/2
c     as, bs, cs, alps, bets, gams - reciprocal lattice parameters

      double precision  a, b, c, alp, bet, gam, v, s, pi, 
     _                  as, bs, cs, alps, bets, gams
      double precision invlen, invang

      pi=3.141593

c           Assign unit cell parameters
c           source:  B. Kahr & R.L. Carter, Mol. Cryst. Liq. Cryst.,
c                       219, 1992, p79.

      a = 14.1562
      b = 21.3190
      c = 13.0654
      alp = 99.9200
      bet = 92.6800
      gam = 106.1500

c           Conversion from degree to radian

      alp = alp * pi / 180.0
      bet = bet * pi / 180.0
      gam = gam * pi / 180.0

c           Calculate reciprocal lattice parameters

      s = (alp+bet+gam)/2.0
      v = 2.0*a*b*c*sqrt(sin(s)*sin(s-alp)*sin(s-bet)*sin(s-gam))
      as = invlen(b, c, alp, v)
      bs = invlen(c, a, bet, v)
      cs = invlen(a, b, gam, v)
      alps = invang(bet, gam, alp)
      bets = invang(gam, alp, bet)
      gams = invang(alp, bet, gam)

      return
      end

      double precision function invlen(x, y, ang, vol)

      double precision x, y, ang, vol

      invlen = x*y*sin(ang)/vol

      return
      end

      double precision function invang(x, y, z)
      double precision x, y, z

      invang = (cos(x)*cos(y)-cos(z))/(sin(x)*sin(y))

      return
      end

      double precision function quad(m, x, y, z, al, be, ga)

      double precision x, y, z, al, be, ga, m(3)
      quad = (m(1)*x)**2+(m(2)*y)**2+(m(3)*z)**2+2.0*m(2)*m(3)* 
     _       y*z*al+2.0*m(3)*m(1)*z*x*be+2.0*m(1)*m(2)*x*y*ga

      return
      end

      double precision function iang(m1,m2,x,y,z,al,be,ga)

      double precision x,y,z,al,be,ga,q1,q2,m1(3),m2(3),pi
      double precision quad

      pi = 3.141593
  
      q1 = quad(m1,x,y,z,al,be,ga)
      q2 = quad(m2,x,y,z,al,be,ga)
      iang = (180.0/pi)*acos((m1(1)*m2(1)*x*x+m1(2)*m2(2)*y*y+
     _       m1(3)*m2(3)*z*z+(m1(2)*m2(3)+m1(3)*m2(2))*y*z*al+
     _       (m1(3)*m2(1)+m1(1)*m2(3))*z*x*be+
     _       (m1(1)*m2(2)+m1(2)*m2(1))*x*y*ga)/sqrt(q1*q2))

      return
      end

      double precision function rmill(a)

      integer a

      rmill=real(a)

      return
      end