Example: Using capping atoms in a periodic system

Here we look at a polyethylene-like (infinite) chain, the carbons being alternatingly substituted with B and N atoms.

/scm-uploads/doc.trunk/Hybrid/_images/BNChain.png

We will use a bunch of different QM systems in a QMMM setup, and check what happens. Inevitably we need to break B-N bonds, and hence capping atoms are used. (In the picture the QM region according to variant one is shown.) We also show the results obtained for the system with the pure QM and MM methods.

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

Download report PeriodicCapping.txt

We optimize the lattice and test several distances

The system can be cut in several variations into a a QM and an MM part breaking a B-N bond

      variation             QM atoms
           var1            B(1),H(2)
           var2            B(5),H(6)
           var3            N(3),H(4)
           var4            N(7),H(8)
           var5  B(1),H(2),N(3),H(4)
           var6  B(5),H(6),N(7),H(8)

Variation one is equivalent to variation two, and variation 3 should be equivalent with variation 4

Variation five is equivalent to variation six

Here are the distances (Angstrom) as obtained with a QM and an MM method
  distance         qm         mm
 B(1)-H(2)      1.182      1.185
 B(5)-H(6)      1.182      1.185
 N(3)-H(4)      1.007      1.045
 N(7)-H(8)      1.007      1.045
 B(1)-N(3)      1.431      1.508
 B(5)-N(7)      1.431      1.508

Now we try the hybrid engine with several variations for the QM region

Two capping methods are tried as well.

      variation         capping          energy       B(1)-H(2)       B(5)-H(6)       N(3)-H(4)       N(7)-H(8)       B(1)-N(3)       B(5)-N(7)
           var1           fixed       -2.901499           1.184           1.185           1.045           1.045           1.508           1.508
           var1      fractional       -2.787165           1.198           1.182           1.044           1.044           1.673           1.505
           var2           fixed       -2.901499           1.185           1.184           1.045           1.045           1.508           1.508
           var2      fractional       -2.787165           1.182           1.198           1.044           1.044           1.505           1.673
           var3           fixed       -4.791110           1.184           1.184           0.993           1.044           1.508           1.507
           var3      fractional       -4.733046           1.184           1.184           0.997           1.045           1.657           1.506
           var4           fixed       -4.791110           1.184           1.184           1.044           0.993           1.507           1.508
           var4      fractional       -4.733046           1.184           1.184           1.045           0.997           1.506           1.657
           var5           fixed       -6.741093           1.189           1.187           1.003           1.045           1.390           1.489
           var5      fractional       -6.648313           1.198           1.183           1.003           1.045           1.405           1.505
           var6           fixed       -6.741093           1.187           1.189           1.045           1.003           1.489           1.390
           var6      fractional       -6.648313           1.183           1.198           1.045           1.003           1.505           1.405

Here are some observations
     * generally the fixed capping seems a bit better

Here are some remarks
     * Starting from the initial very bad structure the fixed capping fails completely for variant 5 and 6
       (not if you use as qm engine band and as mm engine dftb)
     * A reasonable starting geometry can avoid strange collapses
     * The more the two engines disagree about the capped QM region, the stronger the capping forces

Download PeriodicCapping.run

#!/bin/sh

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/var1.xyz
   GuessBonds true
end

Engine DFTB
EndEngine

EOF

aaa1qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#1#2`
aaa2qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#5#6`
bbb1qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#3#4`
bbb2qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#7#8`
ccc1qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#1#3`
ccc2qm=`$AMSBIN/amsreport $AMS_JOBNAME.results/dftb.rkf -r distance#5#7`

printf "We optimize the lattice and test several distances\n" > $report

printf "\nThe system can be cut in several variations into a a QM and an MM part breaking a B-N bond\n" >>$report

printf "\n%15s %20s\n" "variation" "QM atoms" >>$report
printf "%15s %20s\n" "var1" "B(1),H(2)" >>$report
printf "%15s %20s\n" "var2" "B(5),H(6)" >>$report
printf "%15s %20s\n" "var3" "N(3),H(4)" >>$report
printf "%15s %20s\n" "var4" "N(7),H(8)" >>$report
printf "%15s %20s\n" "var5" "B(1),H(2),N(3),H(4)" >>$report
printf "%15s %20s\n" "var6" "B(5),H(6),N(7),H(8)" >>$report

printf "\nVariation one is equivalent to variation two, and variation 3 should be equivalent with variation 4\n" >>$report
printf "\nVariation five is equivalent to variation six\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/var1.xyz
   GuessBonds true
end

Engine ForceField
EndEngine

EOF

aaa1mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#1#2`
aaa2mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#5#6`
bbb1mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#3#4`
bbb2mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#7#8`
ccc1mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#1#3`
ccc2mm=`$AMSBIN/amsreport $AMS_JOBNAME.results/forcefield.rkf -r distance#5#7`

printf "\nHere are the distances (Angstrom) as obtained with a QM and an MM method\n" >> $report
printf "%10s %10s %10s\n"   "distance"  "qm"  "mm" >> $report
printf "%10s %10.3f %10.3f\n"   "B(1)-H(2)"  $aaa1qm $aaa1mm >> $report
printf "%10s %10.3f %10.3f\n"   "B(5)-H(6)"  $aaa2qm $aaa2mm >> $report
printf "%10s %10.3f %10.3f\n"   "N(3)-H(4)"  $bbb1qm $bbb1mm >> $report
printf "%10s %10.3f %10.3f\n"   "N(7)-H(8)"  $bbb2qm $bbb2mm >> $report
printf "%10s %10.3f %10.3f\n"   "B(1)-N(3)"  $ccc1qm $ccc1mm >> $report
printf "%10s %10.3f %10.3f\n"   "B(5)-N(7)"  $ccc2qm $ccc2mm >> $report


printf "\nNow we try the hybrid engine with several variations for the QM region\n" >> $report
printf "\nTwo capping methods are tried as well.\n" >>$report


printf "\n%15s %15s %15s %15s %15s %15s %15s %15s %15s\n" "variation" "capping" "energy" "B(1)-H(2)" "B(5)-H(6)"  "N(3)-H(4)" "N(7)-H(8)" "B(1)-N(3)"  "B(5)-N(7)" >> $report

for system in var1 var2 var3 var4 var5 var6
do

# This calc in only needed to start from a reasonable guess

export AMS_JOBNAME=$system.cheap

rm -rf $AMS_JOBNAME.results

$AMSBIN/ams<<EOF

Task GeometryOptimization

GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=300

System
   GeometryFile $STRUCTDIR/$system.xyz
   GuessBonds true
end

Engine ForceField
EndEngine

EOF


for embedding in 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

GeometryOptimization OptimizeLattice=yes Method=FIRE MaxIterations=100

LoadSystem
 File $system.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
    EndEngine

 EndEngine

EOF

aaa1=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#2`
aaa2=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#5#6`
bbb1=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#3#4`
bbb2=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#7#8`
ccc1=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#1#3`
ccc2=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -r distance#5#7`
xxx=`$AMSBIN/amsreport $AMS_JOBNAME.results/hybrid.rkf -k "AMSResults%Energy"`

printf "%15s %15s %15.6f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f\n" $system $capping $xxx $aaa1 $aaa2 $bbb1 $bbb2 $ccc1  $ccc2 >> $report

done
done
done

printf "\nHere are some observations\n" >>$report
printf "     * generally the fixed capping seems a bit better\n" >>$report

printf "\nHere are some remarks\n" >>$report
printf "     * Starting from the initial very bad structure the fixed capping fails completely for variant 5 and 6\n" >>$report
printf "       (not if you use as qm engine band and as mm engine dftb)\n" >>$report
printf "     * A reasonable starting geometry can avoid strange collapses\n" >>$report
printf "     * The more the two engines disagree about the capped QM region, the stronger the capping forces\n" >>$report

echo "begin report"
cat $report
echo "end report"