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
Task GeometryOptimization
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"