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"