global sys_name_id 
set sys_name_id gly 
set extfold data_bp 

set top ${extfold}/top_all.rtf
set prm ${extfold}/par_all.prm

source ${extfold}/save_${sys_name_id}.chm
set rad [ pb_types2vdwradii atom_types= $types ]

source ${extfold}/qmatoms 
source ${extfold}/outer_mm
source ${extfold}/inner_mm

exec cp ${extfold}/M_rf .
exec cp ${extfold}/charges.c .
############################################################################

  eandg coords= ${sys_name_id}.c  \
	energy= gsbp2.ener \
	gradient= gsbp2.grad \
          theory=hybrid : [ list \
                coupling= shift \
                cutoff= 200 \
                groups = $groups \
                atom_charges= $charges \
                qm_region= $qmatoms \
                gsbp= yes : [ list \
                        restart= yes \
                        inner_region= $inner_mm_numbers \
                        outer_region= $outer_mm_numbers \
                        centre= [ list 1.28123467786000e+00 -2.09759659649645e-01 3.62827519393981e+00 ] \
                        inner_radius= [ expr 9.0 / 0.5292 ] \
                        vdw_radii= $rad \
                        grid_size= 5.0 \
                        inner_grid_size= 2.0 \
                        nr_multipoles= 10 \
                        cir_factor= 2.5 \
                        bc= dhli \
                        spline_order= 3 \
                        eps_mol= 1.0 \
                        eps_solv= 80.0 ] \
                qm_theory= pointcharge : [ list \
                               charges_frag= charges.c \
				 ] \
                mm_theory= dl_poly : [ list \
                               list_option=none \
                               use_pairlist=no \
                               scale14 = { 1.0 1.0 } \
                               atom_types= $types \
                               use_charmm_psf=yes \
                               charmm_psf_file= ${extfold}/${sys_name_id}.psf \
                               charmm_parameter_file= $prm \
                               charmm_mass_file= $top ] ]

set final_energy  [ get_matrix_element matrix=gsbp2.ener indices= {0 0 } ]

# e(mndo99)
set e_ref -2.14157644030298e-01

# Test energy result
if { abs( $final_energy - $e_ref) > 1.0e-7 } {
   chemerr "Failed .. E= $final_energy, should be $e_ref"
}

# Test gradiennt result
subtract_matrix left= gsbp2.grad right= ${extfold}/gsbp2.g.ref result= diff
set max [ get_absmax_matrix_value matrix= diff ]
if { $max > 1.0e-7 } {
   chemerr "Failed .. Maximum absolute deviation of gradient components is $max, should be < 1.0e-7."
} else {
   return OK
}

exec rm M_rf

times
exit