# Example: Using capping atoms in a periodic system with charges¶

Now we look at a BN chain and use charges for the MM calculation.

Let us have a look at the report generated by the example, that pretty much explains what is done

Download report PeriodicCappingWithCharges.txt

We optimize the lattice and test several distances

We divide the system in such a way that there are two equivalent, and hence neutral regions.

Here are the distances (Angstrom) as obtained with a QM and an MM method
distance         qm         mm    err(mm)
B-H      1.182      1.181     -0.001
N-H      1.007      1.041      0.034
B-N      1.431      1.498      0.067

Of course the force field results do not exactly match the QM results, the error displayed in the last column

Now we try the hybrid engine, can we improve the bonds in the QM region?

We start from the geometry calculated with the (cheap) forcefield

In this table we show the errors in bond lengths (in the QM region) of the hybrid method with respect to the QM method

embedding         capping          energy             B-H             N-H             B-N
mechanical           fixed       -6.845015           0.004          -0.007          -0.049
mechanical      fractional       -6.746367           0.013          -0.006          -0.034
electrostatic           fixed       -6.748753           0.002          -0.004          -0.040
electrostatic      fractional       -6.652196           0.010          -0.003          -0.024

Here are some observations
* the B-H distance is a bit worse than with a plain forcefield, especially with fractional capping
* the N-H distance is much better than with the plain forcefield
* the B-N distance is a bit better than with the plain forcefield, now too short. Fractional capping works best.
* Electrostatic embedding is doing slightly better than mechanical embedding, the biggest improvement is on the B-N bond


Download PeriodicCappingWithCharges.run

#!/bin/bash

export NSCM=1

report=report.txt

STRUCTDIR=$AMSHOME/examples/Hybrid/PeriodicCapping/systems # ensure that not a comma is used for decimals in the printf function LC_NUMERIC=en_US.UTF-8 export AMS_JOBNAME=reference rm -rf$AMS_JOBNAME.results

$AMSBIN/ams<<EOF Task GeometryOptimization GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=300 System GeometryFile$STRUCTDIR/var5.xyz
GuessBonds true
end

Engine DFTB
EndEngine

EOF

aaa1qm=$AMSBIN/amsreport$AMS_JOBNAME.results/dftb.rkf -r distance#1#2
bbb1qm=$AMSBIN/amsreport$AMS_JOBNAME.results/dftb.rkf -r distance#3#4
ccc1qm=$AMSBIN/amsreport$AMS_JOBNAME.results/dftb.rkf -r distance#1#3

printf "We optimize the lattice and test several distances\n" > $report printf "\nWe divide the system in such a way that there are two equivalent, and hence neutral regions.\n" >>$report

export AMS_JOBNAME=cheap

rm -rf $AMS_JOBNAME.results$AMSBIN/ams<<EOF

GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=300

System
GeometryFile $STRUCTDIR/var5.xyz LoadForceFieldCharges File=reference.results GuessBonds true end Engine ForceField NonBondedCutoff 50 [Bohr] EndEngine EOF aaa1mm=$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#1#2 bbb1mm=$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#3#4 ccc1mm=$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#1#3 errmma=echo "$aaa1mm- $aaa1qm" | bc errmmb=echo "$bbb1mm- $bbb1qm" | bc errmmc=echo "$ccc1mm- $ccc1qm" | bc printf "\nHere are the distances (Angstrom) as obtained with a QM and an MM method\n" >>$report
printf "%10s %10s %10s %10s\n"   "distance"  "qm"  "mm" "err(mm)">> $report printf "%10s %10.3f %10.3f %10.3f\n" "B-H"$aaa1qm $aaa1mm$errmma >> $report printf "%10s %10.3f %10.3f %10.3f\n" "N-H"$bbb1qm $bbb1mm$errmmb >> $report printf "%10s %10.3f %10.3f %10.3f\n" "B-N"$ccc1qm $ccc1mm$errmmc >> $report printf "\nOf course the force field results do not exactly match the QM results, the error displayed in the last column\n" >>$report

printf "\nNow we try the hybrid engine, can we improve the bonds in the QM region?\n" >> $report printf "\nWe start from the geometry calculated with the (cheap) forcefield\n" >>$report

printf "\nIn this table we show the errors in bond lengths (in the QM region) of the hybrid method with respect to the QM method\n" >> $report printf "\n%15s %15s %15s %15s %15s %15s\n" "embedding" "capping" "energy" "B-H" "N-H" "B-N" >>$report

for system in var5
do
for embedding in  mechanical electrostatic
do

for capping in fixed fractional
do

export AMS_JOBNAME=$system.embedding=$embedding.capping=$capping.go rm -rf$AMS_JOBNAME.results

$AMSBIN/ams<<EOF Task GeometryOptimization Properties Gradients=yes GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=100 LoadSystem File cheap.results End Engine Hybrid Capping AllowHighBondOrders=true Option=$capping

QMMM qmRegion=qm qmEngineID=dftb mmEngineID=ForceField Embedding=$embedding Engine Band EndEngine Engine DFTB EndEngine Engine ForceField NonBondedCutoff 50 [Bohr] EndEngine EndEngine EOF aaa1=$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#2 bbb1=$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#3#4 ccc1=$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#3 xxx=$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k "AMSResults%Energy" # printf "%15s %15s %15s %15.6f %15.3f %15.3f %15.3f\n"$embedding $system$capping $xxx$aaa1 $bbb1$ccc1  >> $report erra=echo "$aaa1- $aaa1qm" | bc errb=echo "$bbb1- $bbb1qm" | bc errc=echo "$ccc1- $ccc1qm" | bc printf "%15s %15s %15.6f %15.3f %15.3f %15.3f\n"$embedding $capping$xxx $erra$errb $errc >>$report

done
done
done

printf "\nHere are some observations\n" >>$report printf " * the B-H distance is a bit worse than with a plain forcefield, especially with fractional capping\n" >>$report
printf "     * the N-H distance is much better than with the plain forcefield \n" >>$report printf " * the B-N distance is a bit better than with the plain forcefield, now too short. Fractional capping works best.\n" >>$report
printf "     * Electrostatic embedding is doing slightly better than mechanical embedding, the biggest improvement is on the B-N bond\n" >>$report echo "begin report" cat$report
echo "end report"