#
#
global root
global theory_type
global qm_theory
global chemsh_default_connectivity_toler 
global chemsh_default_connectivity_scale
global qm
catch {file delete qm_trajectory.xyz}
catch {file delete qm_trajectory.coo}
#
# Future GUI control
#
set groups       undefined
set mxlist       4000
set mxexcl       1000
set diis         1
set scratchdir   .
#
set contyp         undefined
set ctfirst        undefined
set recalc         0
set mm_list_option none
set scfconv        undefined
set ri_memory      1000000
set g98_mem        1000000
set guess          no

# This file was generated by the QUASI - CERIUS2 GUI
# on Thu Feb 15 13:40:33 2001

# Title : Mini TS, QM optimization
set root                   mini_ts
set  conn_scale                1.000000
set  conn_toler                0.500000
set  c2_conn                       0
set theory_type            qm
set function               hdlcopt
set use_zmatrix               0
set residue_treatment          single
set hdlc_region {  } 
set hdlcopt_memory         100
set use_hdlc_constraints   0

# 
# Geometry Optimisation Controls ---------------------------------
# 
set maxstep                200
set newopt_method          bfgs
set find_ts             1
set ts_mode             1
set active_atoms           undefined 

# 
# QM calculation controls ----------------------------------------
# 
set mult                   1
set charge                 -1
set qm_theory              mndo
set basisspec   undefined
set basisfile   undefined
set ecpfile   undefined
set basis   STO3G
set hamiltonian            am1
set scfwf                  rhf
set qm_optstr   undefined

# 
# GAMESS calculation controls ------------------------------------
# 
set direct_scf             1
set symmetry               0
set adaption               0
set use_ri                 0
set qm_restart             0
set maxcyc               200
set qm_accuracy               medium
set chemsh_default_connectivity_toler  0.500000
set chemsh_default_connectivity_scale  1.000000
#
# Whether to pass on the basis set
#
switch $qm_theory {
  mndo - mopac { 
      set bas {}
  }
  default {
      set bas [list basis= $basis ecpfile=$ecpfile  basisfile=$basisfile basisspec = $basisspec ]
  }
}

set newscf {}
# Uncomment here if you need to use the newscf module
# set newscf {  scf_keywords = { { newscf } { maxcyc 100 } {softfail} {end} } }

set qm_theory_args [ list $newscf \
                          hamiltonian = $hamiltonian listing = ${root}_${qm_theory}.log \
	                  accuracy = $qm_accuracy \
                          charge = $charge mult = $mult scftype = $scfwf $bas]
switch $qm_theory {

    gamess {
        set gamess_args [list  direct = $direct_scf symmetry = $symmetry \
		adaption = $adaption maxcyc=$maxcyc ]
	lappend qm_theory_args $gamess_args 
    }
    
    mndo {
	set optarg {}
        if { $qm_optstr != "undefined" } { set optarg "optstr=$qm_optstr" }
        set mndo_args [list idiis=$diis maxcyc=$maxcyc $optarg]
        lappend qm_theory_args $mndo_args 
    }

   turbomole {
        if { $use_ri == "yes" } { 
		set use_ri "use_ri=yes ri_memory=$ri_memory" 
        } else {
		set use_ri "use_ri=no"
	}

	set turbomole_args [list  scratchdir=$scratchdir jobname=$root $use_ri ]
	lappend qm_theory_args $turbomole_args

	# puts stdout [ list $qm_theory_args]
   }

   gaussian {
        # if { $ecpspec != "" } { set ecpspec  "ecpspec=$ecpspec" }
        switch $scfconv {
		undefined { switch $function {
	               		 single_point_energy { set scfmodekey SP }
                       		 default { set scfmodekey Tight }
			    }   
                          }
        	default { set scfmodekey "Conver=$scfconv" }
       }
       if { $scfconv != "undefined" && $maxcyc == 64 } {
          set scfkey "scf=($scfmodekey,maxcycle=$maxcyc)"
       } elseif { $scfconv == "undefined" && $maxcyc != 64 } {
          set scfkey "scf=($scfmodekey,maxcycle=$maxcyc)"
       } elseif { $scfconv != "undefined" && $maxcyc != 64 } {
          set scfkey "scf=($scfmodekey,maxcycle=$maxcyc)"
       } else {
         set scfkey "scf=($scfmodekey,maxcycle=$maxcyc)"
       }
  
#	set gaussian_args [list jobname=$root g98_mem=$g98_mem basis=$basis \
#        $ecpspec guess=$guess $scfkey ]

	set gaussian_args [list jobname=$root g98_mem=$g98_mem basis=$basis \
              guess=$guess $scfkey ]

	lappend qm_theory_args $gaussian_args
   }
}

# these will be overwritten in the hybrid cases
set theory $qm_theory
set theory_args $qm_theory_args
proc save_final_structure {root coords} {
    #
    # Save resulting PDB file
    #
    push_banner_flag 0
    puts stdout "Structure saved as ${root}.pdb"
    write_pdb file=${root}.pdb coords=${coords}

    if { [ llength [ get_cell coords=${coords} ] ] == 9 } {
	puts stdout "Structure saved as ${root}.xtl"
	write_xtl file=${root}.xtl coords=${coords}
    } else {
	puts stdout "Structure saved as ${root}.msi"
	write_msi file=${root}.msi coords=${coords}
    }

    pop_banner_flag 0
}

set act_arg ""
if { "$active_atoms" != "undefined" } { set act_arg "active_atoms = [ list $active_atoms ] "}
proc hdlcopt_update { args } {

    global root
    global theory_type
    global qm_theory

    parsearg update { coords } $args
    write_xyz coords=$coords file=update_${root}.xyz

    switch  $theory_type {
	hybrid {
	    write_xyz coords=hybrid.${qm_theory}.coords file=ttt
	    exec cat ttt >> qm_trajectory.xyz
	    catch {file delete ttt}
	    copy_object type=fragment from=hybrid.${qm_theory}.coords to=ttt
	    exec cat ttt >> qm_trajectory.coo
	    catch {file delete ttt}
	}
    }
    end_module
}

#
# Control parameters
#
#  theory theory_args : 
#  root
#  residue_treatment
#  charmm_pdb_file (only when residue treatment is  )
#

if { $use_hdlc_constraints } { 
    foreach constr $constraints {
	switch [ lindex $constr 0 ] {
	    bond {lappend tt [ lrange $constr 0 2] }
	    angle { lappend tt [ lrange $constr 0 3] }
	    torsion  {lappend tt [ lrange $constr 0 4] }
	}
    }
    set cons_internals "constraints = [ list $tt ] "
} else {
    set cons_internals ""
}


#set reghdl $hdlcopt_memory
set reghdl 0

if { [ get_number_of_atoms coords=${root}.c ] > 1000 } { 
    set cfact 0.5
} else {
    set cfact 0.0
}  

if { "$contyp" == "undefined"  } { set contyp  0 }	
if { "$ctfirst" == "undefined" } { set ctfirst 1 }
if { "$recalc" == "undefined"  } { set recalc  0 }

switch $residue_treatment {

    single {
	set res [ res_selectall coords=${root}.c ]
	set core_atoms [ lindex [ lindex $res 1 ] 0 ]
    }
    pdb {
	set res [ pdb_to_res "$charmm_pdb_file"  ]
	# Residue of interest is assumed to be first
	puts stdout "Warning ... Assume core is 1st residue "
	set core_atoms [ lindex [ lindex $res 0 ] 0 ] 
    }
    select_and_cartesian {
	set res [ list core [list $hdlc_region ] ]
	set core_atoms [ lindex [ lindex $res 0 ] 0 ]
    }
    select_and_pdb {
	#
	#  Combine HDLC core from GUI and 
	#  residues from PDB
	#
	set res [ pdb_to_res "$charmm_pdb_file"  ]

	set r_names [ lindex $res 0 ]
	set r_data  [ lindex $res 1 ]

	catch {unset new_data}
	catch {unset new_names}

	for { set res 0 } { $res < [ llength $r_names ] } { incr res 1 } {

	    set name [ lindex $r_names $res ]
	    set data [ lindex $r_data $res ]
	    
	    set temp {}
	    foreach entry $data {
		if { [ lsearch $hdlc_region $entry ] == -1 } {  lappend temp $entry }
	    }

	    switch [ llength $temp ] {
		0 - 1 - 2 {
		    # the remainder (moving part) is to small to be 
		    # optimised as a DLC residue, by removing the atoms
		    # it will be treated as cartesian 
		}

		default {
		    lappend new_data  $temp
		    lappend new_names $name
		}
	    }
	}

	set res [ list $new_names $new_data ]
	
	set r_core    [ list core [ list  $hdlc_region ] ]
	set res  [ inlist function=merge residues= $res residues2= $r_core ]
	set core_atoms $hdlc_region  
	# puts stdout "Final residue list: $res"
    }
}

set core "core= [list $core_atoms] "
set nvar [ expr 3 * [ llength $core_atoms  ] ]

# TODO - allow mode choice

# puts stdout [ list  lockon=$ts_mode residues= $res  nvar= $nvar $core  \
	recalc=$recalc contyp=$contyp \
	ctfirst=$ctfirst \
	memory=$hdlcopt_memory \
	reghdl=$reghdl \
	cfact=$cfact \
	$act_arg $cons_internals \
	coords=${root}.c \
	result=${root}.ts.optimised \
	maxfun=$maxstep \
	update_procedure=hdlcopt_update \
	theory = $theory : [ list $theory_args ] ]


hdlcopt lockon=$ts_mode residues= $res  nvar= $nvar $core  \
	recalc=$recalc contyp=$contyp \
	ctfirst=$ctfirst \
	memory=$hdlcopt_memory \
	reghdl=$reghdl \
	cfact=$cfact \
	$act_arg $cons_internals \
	coords=${root}.c \
	maxfun=$maxstep \
	result=${root}.ts.optimised \
	update_procedure=hdlcopt_update \
	theory = $theory : [ list $theory_args ]

save_final_structure ${root}_ts ${root}.ts.optimised