C
C $Id: nvtdriver.F,v 1.10 1998/07/16 16:40:08 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      subroutine nvtdriver(ierror)

C     MC simulation

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     locals:

      logical newsub, wrtstuff, ok, rok
      dimension ires(maxres)
      logical picked(maxres), subpicked(maxsub), subpicked1(maxsub)
      character wchar, bchar

CCC---------------------------------------------------------------CCC

C     initialize
      call mcsetup(wchar,bchar,newsub,wrtstuff,ok,rok,
     &     ecrt,dcrt,ires,movedres,imc,igtyp,nstp,iter,ierror)
      if (ierror.ne.0) goto 1000

C     the first MC cycle
      call energy(newsub,eheat,ierror)	
      if(ierror.ne.0) goto 1000
      call atmchg
      call gcartmc(grad)
      if (pme) then
         call pme_directmc(pmedir,pmeself,ecoul,grad)
         if (recipintr) then
            call pme_recip(pmerec,grad)
         else
            call pme_recipmc(pmerec,grad)
         endif
C     .  add PME energies to total energy, subtract classic Coulomb
         elr = pmedir + pmeself + pmerec - ecoul
         etot = etot + elr
C     .  add longe range contribution to heat of formation
         eheat = eheat + 23.061*elr
      endif

      call dovir(vir,grad)
      pressins = denstemp + 0.6947696D5*(vir(4)/boxvol)

      call wrtcoor
      call wrtch
      if (gradient) call wrtgrad(grad,igtyp,gnorm,ierror)
      if (ierror.ne.0) goto 1000
      call printen(imc, movedres, wchar, bchar, ok, rok)
      call mccopynew(ires,movedres,subpicked1)
      naccpt = 1

C     the other MC cycles
      newsub = .false.

      do 200 imc=2,nmc

         call mcupdate(imc,wchar,bchar,npicked,picked,
     &        nsubpicked,subpicked,newsub,wrtstuff,nstp,ierror)

C     .  all energy evaluation are full SCF diagonalizations
C     .  split this up so the convergence crit. can be adapted
C     .  when files (charges!) have to be written

C     .  eecrt = eecrtmc; dencrt = dencrtmc; coordinates/charges are
C     .  NOT written to file

         do 101 j=1,nstp-1

            if (submv) then
               call submove(ires, picked, npicked, movedres,subpicked,
     &              nsubpicked,subpicked1)
               if (movedres.eq.0) goto 101
            else
               call resmove(ires, picked, npicked, movedres,subpicked1)
            endif
            call energy(newsub,eheat,ierror)
            if (ierror.eq.0) then
               if (pme) then
                  call atmchg
                  call gcartmc(grad)
                  call pme_directmc(pmedir,pmeself,ecoul,grad)
                  if (recipintr) then
                     call pme_recip(pmerec,grad)
                  else	
                     call pme_recipmc(pmerec,grad)
                  endif	
C     .           add PME energies to total energy, subtract classic Coulomb
                  elr = pmedir + pmeself + pmerec - ecoul
                  etot = etot + elr
C     .           add longe range contribution to heat of formation
                  eheat = eheat + 23.061*elr
               endif
               call resstep(ires,movedres,ok)
               rok = .true.
            else
               ntry = ntry + 1
               ok = .false.
               rok = .false.
               ierror = 0
            endif
            if (ok) then
               if (.not.pme) then
                  call atmchg
                  call gcartmc(grad)
               endif
               call dovir(vir,grad)
               pressins = denstemp + 0.6947696D5*(vir(4)/boxvol)
               call printen(imc, movedres, wchar, bchar, ok, rok)
               call mccopynew(ires,movedres,subpicked1)
            else

C     .        if configuration is rejected then don't calculate
C     .        the virial in order to save time
               tvir = 0.0
               call printen(imc, movedres, wchar, bchar, ok, rok)
               call mccopyold(ires,movedres)
            endif

            newsub = .false.
 101     continue

C     .  adept convergence criteria
         if (wrtstuff) then
            eecrt = ecrt
            dencrt = dcrt
            wchar = 'W'
         endif

C     .  do the last step

         if (submv) then
            call submove(ires,picked,npicked,movedres,subpicked,
     &           nsubpicked,subpicked1)
            if (movedres.eq.0) goto 110
         else
            call resmove(ires, picked, npicked, movedres,subpicked1)
         endif
         call energy(newsub,eheat,ierror)
         if (ierror.eq.0) then
            if (pme) then
               call atmchg
               call gcartmc(grad)
               call pme_directmc(pmedir,pmeself,ecoul,grad)
               if (recipintr) then
                  call pme_recip(pmerec,grad)
               else
                  call pme_recipmc(pmerec,grad)
               endif
C     .        add PME energies to total energy, subtract classic Coulomb
               elr = pmedir + pmeself + pmerec - ecoul
               etot = etot + elr
C     .        add longe range contribution to heat of formation
               eheat = eheat + 23.061*elr
            endif
            call resstep(ires,movedres,ok)
            rok = .true.
         else
            ntry = ntry + 1
            ok = .false.
            rok = .false.
            ierror = 0
         endif
         if (ok) then
            if (.not.pme) then
               call atmchg
               call gcartmc(grad)
            endif
            call dovir(vir,grad)
            pressins = denstemp + 0.6947696D5*(vir(4)/boxvol)
            call printen(imc, movedres, wchar, bchar, ok, rok)
            call mccopynew(ires,movedres,subpicked1)
         else

C     .     if configuration is rejected then don't calculate
C     .     the virial in order to save time
            tvir = 0.0
            call printen(imc, movedres, wchar, bchar, ok, rok)
            call mccopyold(ires,movedres)
         endif

 110     newsub = .false.

C-RDC C        write coordinates?
         if (wrtstuff) then
            call wrtdiv(imc)
            call wrtcoor
            call wrtch
            if (gradient) call wrtgrad(grad,igtyp,gnorm,ierror)
            if (ierror.ne.0) goto 1000
         endif

C     imc
 200  continue

      call mcclose

 1000 return

      end