Example: Hybrid engine with charged regions
*******************************************

"Download HybridWithCharges.run"

   #!/bin/sh

   # not needed, just slighly faster
   export NSCM=1

   report=report.txt



   printf "Here we treat H3O+ as qm and OH- as the MM region (Optimizing without regions gives two H2O molecules)\n" > $report
   printf "We do this with both mechanical and electrostatic embedding\n" >> $report
   printf  "\n%15s %10s %10s %10s\n" "embedding" "charge" "d(O-O)" "charges" >> $report

   for charge in 0.0 1.0
   do

   if [ "$charge" = "0.0" ]; then
   chargeOinOHm=0.0
   chargeHinOHm=0.0
   chargeOinH3Op=0.0
   chargeHinH3Op=0.0
   else
   chargeOinOHm=-1.123
   chargeHinOHm=0.123
   chargeOinH3Op=-0.5
   chargeHinH3Op=0.5
   fi

   export AMS_JOBNAME=quild.charge=$charge

   rm -rf $AMS_JOBNAME.results

   "$AMSBIN/ams" << eor

   Task GeometryOptimization

   Properties Charges=yes

   GeometryOptimization
     Convergence Gradients=1.0e-6
   End

   System
       Atoms
           O -1.527946410885647 -0.2107366711137158 -0.0008116899510243671 region=QM   ForceField.Charge=$chargeOinH3Op
           H -0.8459142126057956 0.3517312394359257 0.4094504676540848     region=QM   ForceField.Charge=$chargeHinH3Op
           H -1.834953147575289 0.1051014241823828 -0.8704652381864062     region=QM   ForceField.Charge=$chargeHinH3Op
           H -1.328032016244278 -1.164422847242489 0.02894848344144469     region=QM   ForceField.Charge=$chargeHinH3Op
           O 0.6370858511871781 -0.3378071707560572 -0.0006181020627287671 region=MM   ForceField.Charge=$chargeOinOHm
           H 1.318474396634582 0.2241299231185073 0.4092568796869673       region=MM   ForceField.Charge=$chargeHinOHm
       End
       GuessBonds True
   End

   Engine Hybrid
       Energy
          Term Factor=1.0  Region=*  EngineID=ForceField
          Term Factor=-1.0 Region=QM EngineID=ForceField     Charge=$charge
          Term Factor=1.0  Region=QM EngineID=DFTB           Charge=$charge
       End

       Engine DFTB
           Model GFN1-xTB
       EndEngine
       
       Engine ForceField
       EndEngine
   EndEngine

   eor

   ddd=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#5`
   eee=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k AMSResults%Charges#5.3f`

   printf "%15s %10s %10s %10s %10s %10s %10s %10s %10s\n" "mechanical"  $charge $ddd $eee >> $report

   export AMS_JOBNAME=qmmm.charge=$charge

   rm -rf $AMS_JOBNAME.results

   "$AMSBIN/ams" << eor

   Properties Charges=yes

   Task GeometryOptimization

   GeometryOptimization
     Convergence Gradients=1.0e-6
   End

   System
       Atoms
           O -1.527946410885647 -0.2107366711137158 -0.0008116899510243671 region=QM    ForceField.Charge=$chargeOinH3Op
           H -0.8459142126057956 0.3517312394359257 0.4094504676540848     region=QM    ForceField.Charge=$chargeHinH3Op
           H -1.834953147575289 0.1051014241823828 -0.8704652381864062     region=QM    ForceField.Charge=$chargeHinH3Op
           H -1.328032016244278 -1.164422847242489 0.02894848344144469     region=QM    ForceField.Charge=$chargeHinH3Op
           O 0.6370858511871781 -0.3378071707560572 -0.0006181020627287671 region=MM    ForceField.Charge=$chargeOinOHm
           H 1.318474396634582 0.2241299231185073 0.4092568796869673       region=MM    ForceField.Charge=$chargeHinOHm
       End
       GuessBonds True
   End

   Engine Hybrid
       QMMM QMRegion=QM QMEngineID=DFTB MMEngineID=ForceField QMCharge=$charge MMCharge=-$charge

       Engine DFTB
           Model GFN1-xTB
       EndEngine
       
       Engine ForceField
       EndEngine
   EndEngine

   eor

   ddd=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#5`
   eee=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k AMSResults%Charges#5.3f`

   printf "%15s %10s %10s %10s %10s %10s %10s %10s %10s\n" "electrostatic"   $charge $ddd $eee >> $report

   done


   printf "\n* Using charges shortens the O-O distance\n" >> $report
   printf "* In this case the results (mechanical vs. electrostatic) are quite similar as apparently the OH does not polarize the QM region much\n" >> $report



   echo "start of report"
   cat $report
   echo "end of report"






   report=report2.txt

   printf "\nNow we add an extra OH- to the mm region and get a total charge of -1\n" > $report
   printf "We do this with mechanical and electrostatic embedding\n" >> $report
   printf "We look at two distances: d(O1-O5) and d(O1-O7) \n" >> $report
   printf "Atom O1 is in the H3O+ and atoms O5 and O7 are in the two OH- molecules\n" >> $report
   printf  "\n%15s %15s %10s %10s %10s %10s\n" "embedding"  "optim" "d(1-5)" "d(1-7)" "energy" >> $report

   charge=1.0
   chargeOinOHm=-1.123
   chargeHinOHm=0.123
   chargeOinH3Op=-0.5
   chargeHinH3Op=0.5

   for embedding in mechanical electrostatic
   do

   for optim in FIRE # Quasi-Newton
   do

   export AMS_JOBNAME=embedding=$embedding.optim=$optim

   rm -rf $AMS_JOBNAME.results

   "$AMSBIN/ams" << eor

   Task GeometryOptimization

   Properties Charges=yes

   GeometryOptimization
     Method $optim
     MaxIterations 3000
     Convergence Gradients=1.0e-6
   End

   System
       Atoms
           O 0.9019652567984636 -1.133079116834755   0.01338426553857459  region=QM   ForceField.Charge=$chargeOinH3Op
           H 0.1122251167578682 -1.036551903399635   0.5668491423154995   region=QM   ForceField.Charge=$chargeHinH3Op
           H 1.037136681303829  -0.2320347366030556 -0.3773644469587724   region=QM   ForceField.Charge=$chargeHinH3Op
           H 1.678241221654873  -1.266912785246295   0.5779953693196539   region=QM   ForceField.Charge=$chargeHinH3Op
           O -1.130580450693341   0.6009421414132099 -0.02453852439122078  region=MM   ForceField.Charge=$chargeOinOHm
           H -1.671378074377012   1.410809444490273  -0.2141830902463049   region=MM   ForceField.Charge=$chargeHinOHm
           O 3.346891191122751  -0.05485781804516161 0.01059240308504993  region=MM   ForceField.Charge=$chargeOinOHm
           H 4.099773764135065   0.5660034244354222 -0.1683355405307263   region=MM   ForceField.Charge=$chargeHinOHm
       End
       BondOrders
            1 3 1.0
            1 4 1.0
            2 1 1.0
            5 6 1.0
            7 8 1.0
       End
       Charge -1.0
   End

   Engine Hybrid

       QMMM qmRegion=QM qmCharge=1.0 mmCharge=-2.0 qmEngineID=dftb mmEngineID=forcefield Embedding=$embedding

       Engine DFTB
           Model GFN1-xTB
       EndEngine
       
       Engine ForceField
       EndEngine
   EndEngine

   eor


   d15=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#5`
   d17=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#7`
   eee=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k "AMSResults%Energy"`


   printf "%15s %15s %10s %10s %10.4f\n" $embedding $optim $d15 $d17 $eee >> $report

   done
   done


   printf "\n* Very flat PES as function of these two distances\n" >> $report
   printf "\n* Electrostatic embeddiding gives a bit shorter distances\n" >> $report


   echo "start of report"
   cat $report
   echo "end of report"


