# 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

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
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

GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=100

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"
```