! **************************************************************
! This file contains the:  main (SINGLE PROCESSOR JOBS ONLY,
!                                FOR PARALLEL JOBS USE pmain)
!
! Copyright 2003-2005  Frank Eisenmenger, U.H.E. Hansmann,
!                      Shura Hayryan, Chin-Ku
! Copyright 2007       Frank Eisenmenger, U.H.E. Hansmann,
!                      Jan H. Meinke, Sandipan Mohanty
!
! CALLS: init_energy,init_molecule
! CALLS TASK SUBROUTINE: anneal,canon,elp,minim,mulcan_par,
!                        mulcan_sim,partem_s, or regul
! CAN ALSO CALL MEASUREMENT ROUTINES: cnteny,contacts,helix,hbond,
!                                    outpdb,outvar,rgyr,
!                                    rmsinit and rsmdfun,zimmer
! $Id: main.f 334 2007-08-07 09:23:59Z meinke $
! **************************************************************

      program main

      use timer
      include 'INCL.H'
      include 'INCP.H'
      double precision eps, temp, e, energy, enyshe_simple

      integer iabin, imin, maxit, nequi
      integer nsweep, nmes, ncalls, nacalls, i

      common/updstats/ncalls(5),nacalls(5)
      character*80 libdir, seqfile, varfile
      character grpn*4,grpc*4
      logical lrand,bgsposs
      integer argc, status, argv_length
      character(len=255) :: argv

! =================================================== Energy setup

!            Directory for SMMP libraries
!     Change the following directory path to where you want to put SMMP
!     libraries of residues.
      libdir='./SMMP/'

!     Set the maximum log level. The larger the number the more detailed
!     the log.
      MAXLOGLEVEL = 100
!     File unit to use for the log file.
      LOGFILEUNIT = 27
      open(LOGFILEUNIT, file="smmp.log")

!      The switch in the following line is now not used.
      flex=.false.        ! .true. for Flex  / .false. for ECEPP

!     Choose energy type with the following switch instead ...
      ientyp = 0
!        0  => ECEPP2 or ECEPP3 depending on the value of sh2
!        1  => FLEX
!        2  => Lund force field
!        3  => ECEPP with Abagyan corrections
!

      sh2=.false.         ! .true. for ECEPP/2; .false. for ECEPP3
      epsd=.false.        ! .true. for  distance-dependent  dielectric
                          !  permittivity

      itysol= 0    !  0: vacuum
                   ! >0: numerical solvent energy
                   ! <0: analytical solvent energy & gradients

      call init_energy(libdir)

! ================================================= Structure setup

      grpn = 'nh2' ! N-terminal group
      grpc = 'cooh'! C-terminal group

      iabin = 1  ! =0: read from PDB-file
                 ! =1: ab Initio from sequence (& variables)
      seqfile='polyq.seq'
!      seqfile='polyA.pdb'
      varfile='polyq.var'
       varfile = ' '

      ntlml = 0
      write (logString, *) 'Solvent: ', itysol
!     Initialize random number generator.
      call sgrnd(31433)

      if (itysol.eq.0.and.ientyp.eq.3) then
         print *,'Can not use Abagyan entropic corrections without '
         print *,'solvent term. '
         stop
      endif

      call init_molecule(iabin,grpn,grpc,seqfile,varfile)

! Decide if and when to use BGS, and initialize Lund data structures
      bgsprob=0.75   ! Prob for BGS, given that it is possible
! upchswitch= 0 => No BGS 1 => BGS with probability bgsprob
! 2 => temperature dependent choice
      upchswitch=1
      rndord=.true.
      call init_lund
      if (ientyp.eq.2) call init_lundff
      if (ientyp.eq.3) call init_abgn


! ========================================  Add your task down here

      imin = 1 ! Quasi-Newton
      maxit = 15000 ! maximum number of iterations in minimization
      eps = 1.0d-7 ! requested precision
!      call minim(imin, maxit, eps)
!      call outvar(0, ' ')
!     To do a canonical Monte Carlo simulation uncomment the lines below
      nequi = 100
      nsweep = 500
      nmes = 10
      temp = 300.0
      lrand = .true.
      E = energy()
      write (logString, *) E, eyel,eyvw,eyhb,eyvr
      call outpdb(1, "polyrQ.pdb")
!      Canonical Monte Carlo
      call init_timer()
      call start_timer(1)
      call canon(nequi, nsweep, nmes, temp, lrand)
      call stop_timer(1)

      E = 0
!      call start_timer(1)
!      do i =1, nsweep
!         E = E + enyshe_simple()
!      end do
!      call stop_timer(1)
      call evaluate(1)
      print *,timingData(1)%average, "s.", E
!      For simulated annealing uncomment the lines below
!      tmin = 200.0
!      tmax = 500.0
!      call anneal(nequi, nsweep, nmes, tmax, tmin, lrand);
! ========================================  End of main
       end