! DYNAMOCRD_QM_ENERGY_GRADIENT 
! a helper application to get semiempirical QM energy and gradient for a ligand in 
! gelly_refine and autobuster. Uses fdynamo to do this see
! http://www.pdynamo.org/mainpages/downloads.html
!
!===============================================================================
!
!       This program is free software; you can redistribute it and/or
!       modify it under the terms of the GNU General Public License
!       as published by the Free Software Foundation; either version 2
!       of the License, or (at your option) any later version.
!
!       This program is distributed in the hope that it will be useful,
!       but WITHOUT ANY WARRANTY; without even the implied warranty of
!       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!       GNU General Public License for more details.
!
!       You should have received a copy of the GNU General Public License
!       along with this program; if not, write to the Free Software
!       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
!       MA  02110-1301, USA.
!
!===============================================================================
!
! Author: Oliver Smart, Global Phasing Ltd.
! Version: 04/2010
! Email: buster-develop@globalphasing.com
PROGRAM DYNAMOCRD_QM_ENERGY_GRADIENT
USE DYNAMO
IMPLICIT NONE
INTEGER :: IARGC, ARGS_NUMBER, IERR
INTEGER :: MYCHARGE=0, MYMULTIPLICITY=1
CHARACTER(LEN=200)   :: ARGS_STRING, METHOD

WRITE(*,'(a/)')                                                                               &
&'helper: DYNAMOCRD_QM_ENERGY_GRADIENT helper for gelly (used by autobuster refine and grade)'
WRITE(*,'(a/)')                                                                               &
&'PROGRAM_TAG DYNAMOCRD_QM 25.01.11 using FDYNAMOv2.2'

CALL DYNAMO_HEADER
! read in crd - fixed name 
CALL COORDINATES_DEFINE  ( "gelly_coords_for_qmhelper.crdin" )

! get an optional charge and multiplicity on the command line
ARGS_NUMBER = IARGC() ! number of command line arguments
IF (ARGS_NUMBER.GE.1) THEN ! first argument charge
  CALL GETARG( 1, ARGS_STRING)
  READ( ARGS_STRING, *, IOSTAT=IERR) MYCHARGE
  IF (IERR.NE.0) THEN
    WRITE(*,*) 'ERROR reading charge from argument: '//TRIM(ARGS_STRING)
    STOP 'ERROR reading charge'
  ENDIF
ENDIF ! first argument charge

IF (ARGS_NUMBER.GE.2) THEN ! 2nd multipicity
  CALL GETARG( 2, ARGS_STRING)
  READ( ARGS_STRING, *, IOSTAT=IERR) MYMULTIPLICITY
  IF (IERR.NE.0) THEN
    WRITE(*,*) 'ERROR reading multiplicity from argument: '//TRIM(ARGS_STRING)
    STOP 'ERROR reading multiplicity'
  ENDIF
ENDIF 

METHOD="AUTO" ! default
IF (ARGS_NUMBER.GE.3) CALL GETARG( 3, METHOD) ! 3rd arg method
! make case insensitive
IF (METHOD == "am1 ") METHOD = "AM1 "
IF (METHOD == "mndo") METHOD = "MNDO"
IF (METHOD == "pddg") METHOD = "PDDG"
IF (METHOD == "pm3 ") METHOD = "PM3 "
IF (METHOD == "rm1 ") METHOD = "RM1 "
IF (METHOD == "auto") METHOD = "AUTO"
IF (METHOD == "AUTO") THEN
  CALL DYNAMOCRD_QM_AUTO_CHOOSE(METHOD)
ELSEIF ( ( METHOD /= "AM1 " ) .AND. ( METHOD /= "MNDO" ) .AND. &
&        ( METHOD /= "PDDG" ) .AND. ( METHOD /= "PM3 " ) .AND. &
&        ( METHOD /= "RM1 " ) ) THEN
   WRITE(*,*) 'ERROR qm method: '//TRIM(METHOD)//'  not recogized! fdynamo only knows AM1 MNDO PDDG PM3 RM1'
   STOP 'ERROR unrecognized qm method'
ENDIF

WRITE(*,'(a/)')             &
&'METHOD_TAG '//TRIM(METHOD)

WRITE(*,'(/a)')                                                        &
&'helper: increase maximum number of scf cycles allowed from 100 to 300'
CALL MOPAC_SCF_OPTIONS (  ITERATIONS = 300)
WRITE(*,'(/A,I4,A,I4)')                                                                                                 &
&'helper: will now call mopac setup with method='//TRIM(METHOD)//' charge= ', mycharge, ' multiplicity= ', MYMULTIPLICITY
CALL MOPAC_SETUP ( METHOD, CHARGE= MYCHARGE, MULTIPLICITY= MYMULTIPLICITY)
!  find out energy and gradient 
CALL GRADIENT
! write in tagged form
WRITE(*,'(''ENERGYTAG    '',3E30.20)') ETOTAL
WRITE(*,'(''GRADIENTTAG  '',3E30.20)') ATMDER(1:3,1:NATOMS)

END PROGRAM DYNAMOCRD_QM_ENERGY_GRADIENT
!
SUBROUTINE DYNAMOCRD_QM_AUTO_CHOOSE(METHOD)
USE DYNAMO
IMPLICIT NONE
CHARACTER(LEN=200), INTENT(OUT) ::  METHOD
! nasty hack to choose best qm method that covers all elements 
! got list of which elements covered by each method
! egrep "SEPAR|PARAMETERS_" fdynamo-2.2/source/mopac_parameters.F90 > temp.F90
LOGICAL, DIMENSION(1:NELEMENTS) :: mySEPAR

METHOD='NONE'
IF (METHOD=='NONE') THEN ! RM1?
      mySEPAR(:) = .FALSE.
      mySEPAR ( 1) = .TRUE.
      mySEPAR ( 6) = .TRUE.
      mySEPAR ( 7) = .TRUE.
      mySEPAR ( 8) = .TRUE.
      mySEPAR ( 9) = .TRUE.
      mySEPAR (15) = .TRUE.
      mySEPAR (16) = .TRUE.
      mySEPAR (17) = .TRUE.
      mySEPAR (35) = .TRUE.
      mySEPAR (53) = .TRUE.
      IF (ALL(mySEPAR(ATMNUM(1:NATOMS)))) METHOD='RM1'
ENDIF

IF (METHOD=='NONE') THEN ! PM3?
      mySEPAR(:) = .FALSE.
      mySEPAR(1)   = .TRUE.
      mySEPAR(4)   = .TRUE.
      mySEPAR(6)   = .TRUE.
      mySEPAR(7)  = .TRUE.
      mySEPAR(8)   = .TRUE.
      mySEPAR(9)   = .TRUE.
      mySEPAR(11) = .TRUE.
      mySEPAR(12) = .TRUE.
      mySEPAR(13)  = .TRUE.
      mySEPAR(14)  = .TRUE.
      mySEPAR(15)  = .TRUE.
      mySEPAR(16)  = .TRUE.
      mySEPAR(17)  = .TRUE.
      mySEPAR(19) = .TRUE.
      mySEPAR(30)  = .TRUE.
      mySEPAR(31)  = .TRUE.
      mySEPAR(32) = .TRUE.
      mySEPAR(33)  = .TRUE.
      mySEPAR(34)  = .TRUE.
      mySEPAR(35)  = .TRUE.
      mySEPAR(48)  = .TRUE.
      mySEPAR(49)  = .TRUE.
      mySEPAR(50) = .TRUE.
      mySEPAR(51) = .TRUE.
      mySEPAR(52) = .TRUE.
      mySEPAR(53)  = .TRUE.
      mySEPAR(80)  = .TRUE.
      mySEPAR(81)  = .TRUE.
      mySEPAR(82)  = .TRUE.
      mySEPAR( 83) = .TRUE.
      IF (ALL(mySEPAR(ATMNUM(1:NATOMS)))) METHOD='PM3'
ENDIF


IF (METHOD=='NONE') THEN ! AM1?
      mySEPAR(1) = .TRUE.
      mySEPAR(3) = .TRUE.
      mySEPAR(4) = .TRUE.
      mySEPAR(5) = .TRUE.
      mySEPAR(6) = .TRUE.
      mySEPAR(7) = .TRUE.
      mySEPAR(8) = .TRUE.
      mySEPAR(9) = .TRUE.
      mySEPAR(11) = .TRUE.
      mySEPAR(12) = .TRUE.
      mySEPAR(13) = .TRUE.
      mySEPAR(14) = .TRUE.
      mySEPAR(15) = .TRUE.
      mySEPAR(16) = .TRUE.
      mySEPAR(17) = .TRUE.
      mySEPAR(19) = .TRUE.
      mySEPAR(30)  = .TRUE.
      mySEPAR(35) = .TRUE.
      mySEPAR(53) = .TRUE.
      IF (ALL(mySEPAR(ATMNUM(1:NATOMS)))) METHOD='AM1'
ENDIF

IF (METHOD=='NONE') THEN 
   WRITE(*,'(A)')                                                            &
   &'helper: ERROR no QM method in fdynamo cannot cope with all your elements!'
   STOP 'ERROR no QM method in fdynamo cannot cope with all your elements!'
ENDIF


WRITE(*,'(A)')                                    &
&'helper: AUTO option choose method '//TRIM(METHOD)
RETURN
END SUBROUTINE DYNAMOCRD_QM_AUTO_CHOOSE
