#!/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"