      program nmr_fit
C     ===============

C--------------------------------------------------------------------------
C
C     Program designed to use a Monte Carlo method to fit Si-NMR spectra
C
C     Version 1.1 (first working version)
C     Martin Dove, February 1995
C
C--------------------------------------------------------------------------
C
C
C     This program is available for free use. It can also be freely
C     distributed, although it would be useful if any new users could
C     e-mail the author (martin @ minp.esc.cam.ac.uk) to be added to a
C     users' list, not least in case there is a need to provide information
C     about upgrades or (if we are unlucky) any mistakes that have been
C     fixed.
C
C     Please send any suggestions, comments or corrections by e-mail to
C     martin @ minp.esc.cam.ac.uk
C
C     Any acknowledgement of the use of this program should be made by
C     reference to the paper that describes the program:
C
C         "The use of Monte Carlo methods to determine the distribution of
C         Al and Si cations in framework aluminosilicates from 29-Si
C         MAS-NMR data"
C
C         M T Dove and V Heine
C
C         American Mineralogist (submitted)
C
C
C--------------------------------------------------------------------------
C
C                           INFORMATION
C                           ===========
C
C     Store of information
C     ~~~~~~~~~~~~~~~~~~~~
C
C     Information about each atom is stored in the character string
C     cneigh(natoms). This string contains 6 integers. Integer 1 is the
C     atom type (0 = Si, 1 = Al), Integer 2 gives the type of site, and
C     Integers 3-6 give the numbers of the connected atoms.
C
C
C     Input files
C     ~~~~~~~~~~~
C
C     Information is read from 3 files:
C
C     Channel 11 has the information about the crystal structure:
C           Line  1  = the space group symbol
C           Line  2  = the six cell parameters in the second line
C           Lines 3+ = a number =1 for Al and =0 for Si, a number to label
C                      the site and the fractional coordinates - one 
C                      line per symmetrically-independent atom. The first
C                      parameter is only relevant if the number of
C                      randomisation steps (channel 12) is non-zero. If it
C                      is zero you are starting from an ordered structure,
C                      otherwise you are starting from random order.
C
C     Channel 12 contains information about the simulation:
C           Line  1  = the number of unit cells required along the x, y and 
C                      z directions
C           Line  2  = the number of initial randomisation steps
C           Line  3  = the relative proportions of Al and Si atoms
C           Lines 4+ = the number of Monte Carlo steps required, the 
C                      parameter beta which as as the reciprocal of a given
C                      "temperature", and a number to indicate the average
C                      number of steps between updates of the accumulated
C                      results. This last line will be repeated until the
C                      end of the file, to allow you to "cool down".
C
C     Channel 13 contains information about the nmr data:
C           Lines 1+ = the first parameter is the label of the site, then
C                      this is followed by the number of sites with 0, 1, 2
C                      3 and 4 Al neighbours. The final parameter (#6) is
C                      an error for the nmr data, assumed to be the same 
C                      for each of the nmr peaks.
C
C
C     Ouput files
C     ~~~~~~~~~~~
C
C     There are a number of output files, which are not named but simply
C     written to files of the form (at least in VMS) FORxxx.DAT:
C
C     Channel 8:   Information about the generation of the structure
C     Channel 9:   List of atoms (fractional and orthorgonal units)
C     Channel 21+: List of numbers of Al-Al bonds for each site, one line
C                  per MC analysis. Each file is for a different type of
C                  site
C     Channel 41+: List of nmr peak heights (normalised so that the sum is
C                  unity). Each line has a value of Si-nAl for each value 
C                  of n (0-4), one line per MC analysis. Each file is for a
C                  different type site.
C     Channel 61:  List of number of Al per site, one line per MC analysis.
C     Channel 71:  Agreement factor for each MC analysis.
C     Channel 81:  List of all sites with atom types and neighbours
C
C
C     Some programming details
C     ~~~~~~~~~~~~~~~~~~~~~~~~
C
C     The program uses three NAG routines for the random number generators.
C
C     integer function g05dyf(m,n): gives a random integer between m and
C                                   n selected from a uniform distribution
C
C     subroutine g05ccf:            selects a random seed for the random
C                                   number generators
C
C     real*8 function g05caf(x):    gives a random number (real*8) selected
C                                   from a uniform distribution between
C                                   0.0 and 1.0. x is a dummy real*8
C                                   variable.
C
C     If NAG is not available, these can be substituted by standard
C     alternatives. Indeed, alternatives are available in an accompanying
C     file.
C
C--------------------------------------------------------------------------

      common /symmetry/ lattice_info(2),lattice_type,operators(96)
      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /numbers/ nal(12),nalal(12,12),nsinal(12,5)

      character*80 operators,cneigh,cbuffer
      character*1 lattice_type
      character spg(20)*1,latsym*1

      read(11,'(a)') cbuffer

      if (cbuffer(1:7).eq.'restart') then
        call restart
        go to 1001
      endif

      read(cbuffer,'(20a1)') spg

C ... call spacegroup program to generate all space group operations
      call spaceg(spg)
      write(*,*) ' Space group symmetry generated'

C ... call subroutine to generate structure, reading structure from file
C     and constructing all fractional and orthogonal coordinates
      call generate_structure
      write(*,*) ' Structure generated'

C ... call subroutine to create supercell
      call expand_cell
      write(*,*) ' Structure expanded to form supercell'

C ... call subroutine to establish all neighbours for each T site
      call setup_neighbours
      write(*,*) ' Neighbour lists constructed'

C ... call subroutine to randomise the Al over the T sites
      call randomize
      write(*,*) ' Random assignment of initial Al sites completed'

C ... call subroutine to perform initial analysis of configuration
 1001 call initial_analysis
      write(*,*) ' Initial analysis performed'

      write(*,*) ' Monte Carlo steps now beginning...'

C ... call subroutine to do the Monte Carlo stuff
      call monte_carlo

C ... Finally, write a restart file
      write(81,'(a)') 'restart'
      do 1002 i = 1,natoms
 1002 write(81,'(a)') cneigh(i)

      end


C==========================================================================


      subroutine restart
C     ==================

C--------------------------------------------------------------------------
C
C     Restarts from a previous configuration
C
C--------------------------------------------------------------------------

      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)

      character*80 cneigh,cbuffer

      natoms = 0
 1001 read(11,'(a)',end=2000) cbuffer
      natoms = natoms + 1
      cneigh(natoms) = cbuffer
      go to 1001

 2000 continue

      return
      end


C==========================================================================


      subroutine generate_structure
C     =============================

C--------------------------------------------------------------------------
C
C     This subroutine generates the crystal structure from the space
C     group operations and the input atomic coordinates.
C
C--------------------------------------------------------------------------


      common /lattice_parameters/ a,b,c,al,be,ga
      common /lattice_vectors/ cell(3,3)
      common /symmetry/ lattice_info(2),lattice_type,operators(96)
      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)

      character*80 operators,cneigh
      character*1 lattice_type
      character*2 cx1,cx2,cy1,cy2,cz1,cz2
      logical okay

      data cell/9*0.0/

C ... first form the lattice vectors
      read(11,*) a,b,c,al,be,ga
      al=al/57.295779512
      be=be/57.295779512
      ga=ga/57.295779512
      s=0.5*(al+be+ga)
      vol=2.0*a*b*c*sqrt(sin(s)*sin(s-al)*sin(s-be)*sin(s-ga))
      cell(1,1)=vol/(sin(al)*b*c)
      cell(2,1)=a*cos(ga)*sin(al)
      cell(3,1)=a*cos(be)
      cell(1,2)=0.0
      cell(2,2)=b*sin(al)
      cell(3,2)=b*cos(al)
      cell(1,3)=0.0
      cell(2,3)=0.0
      cell(3,3)=c

C ... read in the fractional coordinates, generate all the related 
C     coordinates by looping over symmetry operations, the centre, and the
C     lattice types, and then check that there are no duplicates.

      natoms = 0
 5000 read(11,*,end=1000) jal,jtype,xr,yr,zr
      if (jal.ne.1) jal = 0
      do 1010 i = 1,lattice_info(1)
      read(operators(i),'(5x,3(i1,1x,i1,2a2,2x))')ix1,ix2,cx1,cx2,
     +            iy1,iy2,cy1,cy2,iz1,iz2,cz1,cz2 
      x = 0.0
      y = 0.0
      z = 0.0
      if (ix2.ne.0) x = float(ix1)/float(ix2)
      if (iy2.ne.0) y = float(iy1)/float(iy2)
      if (iz2.ne.0) z = float(iz1)/float(iz2)
      if (cx1.eq.'+x') x = x + xr
      if (cx1.eq.'-x') x = x - xr
      if (cx1.eq.'+y') x = x + yr
      if (cx1.eq.'-y') x = x - yr
      if (cx1.eq.'+z') x = x + zr
      if (cx1.eq.'-z') x = x - zr
      if (cy1.eq.'+x') y = y + xr
      if (cy1.eq.'-x') y = y - xr
      if (cy1.eq.'+y') y = y + yr
      if (cy1.eq.'-y') y = y - yr
      if (cy1.eq.'+z') y = y + zr
      if (cy1.eq.'-z') y = y - zr
      if (cz1.eq.'+x') z = z + xr
      if (cz1.eq.'-x') z = z - xr
      if (cz1.eq.'+y') z = z + yr
      if (cz1.eq.'-y') z = z - yr
      if (cz1.eq.'+z') z = z + zr
      if (cz1.eq.'-z') z = z - zr
      if (cx2.eq.'+x') x = x + xr
      if (cx2.eq.'-x') x = x - xr
      if (cx2.eq.'+y') x = x + yr
      if (cx2.eq.'-y') x = x - yr
      if (cx2.eq.'+z') x = x + zr
      if (cx2.eq.'-z') x = x - zr
      if (cy2.eq.'+x') y = y + xr
      if (cy2.eq.'-x') y = y - xr
      if (cy2.eq.'+y') y = y + yr
      if (cy2.eq.'-y') y = y - yr
      if (cy2.eq.'+z') y = y + zr
      if (cy2.eq.'-z') y = y - zr
      if (cz2.eq.'+x') z = z + xr
      if (cz2.eq.'-x') z = z - xr
      if (cz2.eq.'+y') z = z + yr
      if (cz2.eq.'-y') z = z - yr
      if (cz2.eq.'+z') z = z + zr
      if (cz2.eq.'-z') z = z - zr
        do 1001 j = 1,(lattice_info(2)+1)
          if (j.eq.2) then
            x = -x
            y = -y
            z = -z
          endif
          if (x.lt.0.0) x = x + 1.0
          if (y.lt.0.0) y = y + 1.0
          if (z.lt.0.0) z = z + 1.0
          if (x.gt.1.0) x = x - 1.0
          if (y.gt.1.0) y = y - 1.0
          if (z.gt.1.0) z = z - 1.0
          lattice = 2
          if (lattice_type.eq.'P') lattice = 1
          if (lattice_type.eq.'F') lattice = 4
          do 1020 k = 1,lattice
            xff = x
            yff = y
            zff = z
            if ((k.eq.2).and.(lattice_type.eq.'A')) then
              yff = y + 0.5
              zff = z + 0.5
            endif
            if ((k.eq.2).and.(lattice_type.eq.'B')) then
              xff = x + 0.5
              zff = z + 0.5
            endif
            if ((k.eq.2).and.(lattice_type.eq.'C')) then
              xff = x + 0.5
              yff = y + 0.5
            endif
            if ((k.eq.2).and.(lattice_type.eq.'I')) then
              xff = x + 0.5
              yff = y + 0.5
              zff = z + 0.5
            endif
            if ((k.eq.2).and.(lattice_type.eq.'F')) then
              xff = x + 0.5
              yff = y + 0.5
            endif
            if ((k.eq.3).and.(lattice_type.eq.'F')) then
              xff = x + 0.5
              zff = z + 0.5
            endif
            if ((k.eq.4).and.(lattice_type.eq.'F')) then
              yff = y + 0.5
              zff = z + 0.5
            endif
            if (xff.ge.1.0) xff = xff - 1.0
            if (yff.ge.1.0) yff = yff - 1.0
            if (zff.ge.1.0) zff = zff - 1.0
            okay = .true.
            if (natoms.gt.0) then
              do 1030 ij = 1,natoms
              if ( (abs(xff-xf(ij)).lt.0.001)
     +          .and. (abs(yff-yf(ij)).lt.0.001)
     +          .and. (abs(zff-zf(ij)).lt.0.001) ) okay = .false.
 1030         continue
            endif
            if (okay) then
              natoms = natoms + 1
              xf(natoms) = xff
              yf(natoms) = yff
              zf(natoms) = zff
              xo(natoms) = cell(1,1)*xff + cell(1,2)*yff + cell(1,3)*zff
              yo(natoms) = cell(2,1)*xff + cell(2,2)*yff + cell(2,3)*zff
              zo(natoms) = cell(3,1)*xff + cell(3,2)*yff + cell(3,3)*zff
              write(cneigh(natoms),*) jal,jtype,0,0,0,0
            endif
 1020     continue
 1001   continue
 1010 continue
      go to 5000
 1000 continue

      write(8,*) ' Number of sites: ',natoms
      return
      end


C==========================================================================


      subroutine expand_cell
C     ======================

C--------------------------------------------------------------------------
C
C     This subroutine expands the list of atoms by adding extra lattice
C     vectors to the basis set.
C
C--------------------------------------------------------------------------

      common /lattice_vectors/ cell(3,3)
      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /expand/ ncell(3)

      character*80 cneigh
      dimension info(6)

      do 1101 i = 1,natoms
 1101 write(9,'(i5,6f10.5)') i,xf(i),yf(i),zf(i),xo(i),yo(i),zo(i)

      read(12,*) ncell 
      n = natoms
      do 1001 i = 1,ncell(1)
      do 1001 j = 1,ncell(2)
      do 1001 k = 1,ncell(3)
      if ((i.eq.1).and.(j.eq.1).and.(k.eq.1)) go to 1001
      ii = i-1
      jj = j-1
      kk = k-1
      do 1002 ia = 1,natoms
      n = n + 1
      xf(n) = xf(ia) + ii
      yf(n) = yf(ia) + jj
      zf(n) = zf(ia) + kk
      xo(n) = xo(ia) + ii*cell(1,1) + jj*cell(1,2) + kk*cell(1,3)
      yo(n) = yo(ia) + ii*cell(2,1) + jj*cell(2,2) + kk*cell(2,3)
      zo(n) = zo(ia) + ii*cell(3,1) + jj*cell(3,2) + kk*cell(3,3)
      write(9,'(i5,6f10.5)') n,xf(n),yf(n),zf(n),xo(n),yo(n),zo(n)
      cneigh(n) = cneigh(ia)
 1002 continue
 1001 continue
      do 1003 i = 1,3
      do 1003 j = 1,3
 1003 cell(i,j) = cell(i,j)*ncell(j)
      natoms = natoms*ncell(1)*ncell(2)*ncell(3)
      if (n.ne.natoms) then
        write(6,*) ' **** Something is wrong with the expansion of the',
     +             ' cell'
        stop 'Marker C'
      endif
      return
      end      


C==========================================================================


      subroutine setup_neighbours
C     ===========================

C--------------------------------------------------------------------------
C
C     This subroutine sorts out the neighbour lists.
C
C--------------------------------------------------------------------------

      common /lattice_parameters/ a,b,c,al,be,ga
      common /lattice_vectors/ cell(3,3)
      common /symmetry/ lattice_info(2),lattice_type,operators(96)
      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /expand/ ncell(3)

      dimension info(6)
      character*80 cneigh

      do 2001 i = 1,natoms
      in = 0
      read(cneigh(i),*) info
      do 2000 j = 1,natoms
      rx = xo(i) - xo(j)
      ry = yo(i) - yo(j)
      rz = zo(i) - zo(j)
      do 1008 l = -1,1
      do 1008 m = -1,1
      do 1008 n = -1,1
      rxx = rx + l*cell(1,1) + m*cell(1,2) + n*cell(1,3)
      ryy = ry + l*cell(2,1) + m*cell(2,2) + n*cell(2,3)
      rzz = rz + l*cell(3,1) + m*cell(3,2) + n*cell(3,3)
      r = sqrt(rxx**2 + ryy**2 + rzz**2)
      if ((r.lt.3.4).and.(r.gt.0.1)) then
        in = in + 1
        if (in.gt.4) then
           write(6,*) '*** Number of neighbours is too great!'
           stop 'Marker D'
        endif
        info(in+2) = j
      endif
 1008 continue
 2000 continue
 2001 write(cneigh(i),*) info

      write(9,*)
      do 2101 i = 1,natoms
 2101 write(9,'(a)') cneigh(i)
      return
      end


C==========================================================================


      subroutine randomize
C     ====================

C--------------------------------------------------------------------------
C
C     This subroutine decides how many Al to use, and then sets the first
C     block of atoms to be Al. Then it chooses pairs of atoms randomly
C     and simply swaps them around.
C
C     Finally, it analyses the final configuration. It calculates the 
C     following sums:
C
C          nal(i)      =  the number of Al on site of type i
C          nalal(i,j)  =  the number of Al-Al linkages between sites i and j
C          nsinal(i,n) =  the number of Si on site type i that are 
C                         linked to n Al
C
C--------------------------------------------------------------------------

      common /lattice_parameters/ a,b,c,al,be,ga
      common /lattice_vectors/ cell(3,3)
      common /symmetry/ lattice_info(2),lattice_type,operators(96)
      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /numbers/ nal(12),nalal(12,12),nsinal(12,5)

      integer g05dyf
      dimension info(6),jnfo(6)
      character*80 cneigh

C ... read in the number of randomisation steps
      read(12,*) nrand

C ... Leave if you are aiming to start from an ordered structure.
      if (nrand.eq.0) then
        return
      endif

C ... read in the relative proportions of Al to Si
      read(12,*) fal, fsi
      nals = nint(natoms*fal/(fal+fsi))

      do 1001 i = 1,nals
      read(cneigh(i),*) info
      info(1) = 1
 1001 write(cneigh(i),*) info

C ... set up a random seed for the random number generator
      call g05ccf

C ... Now mix up the atoms
      do 1002 irand = 1,nrand
       iatom = g05dyf(1,natoms)
       if ((iatom.gt.natoms).or.(iatom.lt.1)) then
        write(6,*) '**** An atom has been selected that does not exist'
        stop 'Marker A'
       endif
 1101  jatom = g05dyf(1,natoms)
       if ((jatom.gt.natoms).or.(jatom.lt.1)) then
        write(6,*) '**** An atom has been selected that does not exist'
        stop 'Marker B'
       endif
       read(cneigh(iatom),*) info
       read(cneigh(jatom),*) jnfo
       if (info(1).eq.jnfo(1)) go to 1101
       jj = jnfo(1)
       jnfo(1) = info(1)
       info(1) = jj
       write(cneigh(iatom),*) info
       write(cneigh(jatom),*) jnfo
 1002 continue

C ... Check that we have a properly averaged number of Al in the first
C     NALS sites
      n = 0
      do 1102 i = 1,nals
      read(cneigh(i),*) j
 1102 n = n + j
      write(8,*) ' Number of Al in the corresponding number of site: ',n

      return
      end


C==========================================================================


      subroutine initial_analysis
C     ===========================

C--------------------------------------------------------------------------
C
C     we need to perform an initial analysis for the initial MC step
C
C--------------------------------------------------------------------------

      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /numbers/ nal(12),nalal(12,12),nsinal(12,5)

      dimension info(6),jnfo(6)
      character*80 cneigh

C ... Now analyse configuration produced, starting with initialising of
C     arrays used to store numbers
      do 1004 i = 1,12
       nal(i) = 0
      do 1005 j = 1,5
 1005 nsinal(i,j) = 0
      do 1004 j = 1,12
       nalal(i,j) = 0
 1004 continue

C ... First, the number of Al on each type of site
      do 2000 i = 1,natoms
       read(cneigh(i),*) info
       nal(info(2)) = nal(info(2)) + info(1)
 2000 continue

C ... Second, the number of Al...Al contacts
      do 2001 i = 1,natoms
       read(cneigh(i),*) info
       do 2001 j = 3,6
        read(cneigh(info(j)),*) jnfo
        nalal(info(2),jnfo(2)) = nalal(info(2),jnfo(2)) + 
     +                           info(1)*jnfo(1)
        nalal(jnfo(2),info(2)) = nalal(jnfo(2),info(2)) + 
     +                           info(1)*jnfo(1)
 2001 continue

C ... Third, evaluate the number of Al atoms about each Si
      do 2002 i = 1,natoms
       read(cneigh(i),*) info
       if (info(1).eq.1) go to 2002
       m = 1
       do 2003 j = 3,6
       read(cneigh(info(j)),*) ja
 2003  m = m + ja
       if (m.gt.5) then
         write(6,*) ' **** Warning, too many Al about a Si'
         stop 'Marker E'
       endif
       nsinal(info(2),m) = nsinal(info(2),m) + 1
 2002 continue

      return
      end


C==========================================================================


      subroutine monte_carlo
C     ======================

C--------------------------------------------------------------------------
C
C     This subroutine does the Monte Carlo rearrangements
C
C--------------------------------------------------------------------------

      common /coordinates/ xf(2000),yf(2000),zf(2000),xo(2000),yo(2000),
     +                     zo(2000),natoms
      common /neighbours/ cneigh(2000)
      common /numbers/ nal(12),nalal(12,12),nsinal(12,5)

      integer g05dyf
      double precision g05caf
      double precision dummy
      dimension info(6),jnfo(6),msinal(12,5),gsinal(12,5)
      character*80 cneigh,cbuffer
      logical okay,nmrdata(12)

C ... set up a random seed for the random number generator
      call g05ccf

      als = 0.0
      do 1998 i = 1,12
      als = als + nal(i)
 1998 nmrdata(i) = .false.
      sils = natoms - als

C ... read in the experimental data
 1999 read(13,'(a)',end=1000) cbuffer
      read(cbuffer,*) i
      nmrdata(i) = .true.
      read(cbuffer,*) junk,(gsinal(i,j),j=1,5),sigma
      go to 1999

C ... read in the number of Monte_Carlo steps and inverse temperature
 1000 read(12,*,end=9000) nmcs,beta,iconfig
      ianal = 0

C ... Loop over the Monte Carlo steps
      do 2000 imcs = 1,nmcs

 1102  iatom = g05dyf(1,natoms)
       if ((iatom.gt.natoms).or.(iatom.lt.1)) then
        write(6,*) '**** An atom has been selected that does not exist'
        stop 'Marker H'
       endif
 1101  jatom = g05dyf(1,natoms)
       if ((jatom.gt.natoms).or.(jatom.lt.1)) then
        write(6,*) '**** An atom has been selected that does not exist'
        stop 'Marker I'
       endif
       read(cneigh(iatom),*) info
       read(cneigh(jatom),*) jnfo
       if (info(1).eq.jnfo(1)) go to 1101

C ... What I do now is the most inefficient thing possible. But, it is
C     probably better to waste the computer's time than my time!       
C     First switch the atoms in any case
       jj = jnfo(1)
       jnfo(1) = info(1)
       info(1) = jj
       write(cneigh(iatom),*) info
       write(cneigh(jatom),*) jnfo

C ... Then obtain the new distribution
       do 1002 i = 1,12
       do 1002 j = 1,5
 1002  msinal(i,j) = 0
       do 2002 i = 1,natoms
        read(cneigh(i),*) info
        if (info(1).eq.1) go to 2002
        m = 1
        do 2003 j = 3,6
         read(cneigh(info(j)),*) ja
 2003   m = m + ja
        if (m.gt.5) then
          write(6,*) ' **** Warning, too many Al about a Si'
          stop 'Marker J'
        endif
        msinal(info(2),m) = msinal(info(2),m) + 1
 2002  continue

C ... Calculate the new and old "energies"
      energyold = 0.0
      energynew = 0.0
      do 2008 i = 1,12
      if (nmrdata(i)) then
        do 2007 j = 1,5
         gn = float(nsinal(i,j))/sils
         gm = float(msinal(i,j))/sils
         energyold = energyold + (gsinal(i,j)-gn)**2
         energynew = energynew + (gsinal(i,j)-gm)**2
 2007   continue
      endif
 2008 continue
      if (energynew.le.energyold) then
        okay = .true.
      else
        random = g05caf(dummy)
        e = (energynew - energyold)/sigma**2
        test = exp(-beta*e)
        if (random.le.test) then
           okay = .true.
        else
           okay = .false.
        endif
      endif

C ... If we like the flip, update the neighbour lists and carry on
      if (okay) then
        do 2004 i = 1,12
        do 2004 j = 1,5
 2004   nsinal(i,j) = msinal(i,j)
        go to 2100
      endif

C ... So we are only here if we don't like the flip, so we now have to
C     set everything back to normal. Afterwards, don't count this is as a
C     step but go back to the random number generators

       read(cneigh(iatom),*) info
       read(cneigh(jatom),*) jnfo
       jj = jnfo(1)
       jnfo(1) = info(1)
       info(1) = jj
       write(cneigh(iatom),*) info
       write(cneigh(jatom),*) jnfo
       go to 1102

C ... We arrive here if we have updated the configuration. We now want to
C     add the results to the store of the results.

 2100  ianal = ianal + 1
       if (ianal.ne.iconfig) go to 2000

       ianal = 0
       write(71,*) energynew

C ... First, the number of Al on each type of site
       do 3101 i = 1,12
 3101  nal(i) = 0
       do 3100 i = 1,natoms
        read(cneigh(i),*) info
        nal(info(2)) = nal(info(2)) + info(1)
 3100  continue
       write(61,'(12i6)') nal

C ... Second, the number of Al-Al contacts
       do 3000 i = 1,12
       do 3000 j = 1,12
 3000  nalal(i,j) = 0
       do 3001 i = 1,natoms
        read(cneigh(i),*) info
        do 3001 j = 3,6
         read(cneigh(info(j)),*) jnfo
         nalal(info(2),jnfo(2)) = nalal(info(2),jnfo(2)) + 
     +                            info(1)*jnfo(1)
         nalal(jnfo(2),info(2)) = nalal(jnfo(2),info(2)) + 
     +                            info(1)*jnfo(1)
 3001  continue

C ... Third, the calculated nmr spectra (i.e. relative heights of peaks
C     from Si with N Al neighbours)
       do 3002 i = 1,12
       if (nmrdata(i)) then
         iout = i + 20
         write(iout,'(12i6)') (nalal(i,j),j=1,12)
         iout = i + 40
         write(iout,'(5f10.5)') ((float(nsinal(i,j))/sils),j=1,5)
       endif
 3002  continue

 2000 continue
      write(*,*) ' Block of Monte Carlo steps completed'
      go to 1000

 9000 continue
      return
      end
