Example: Molecular gun with the hybrid engine¶
In this example we are going to really stretch the use of the Hybrid Engine, and shoot bullets (treated with a QM engine) at a surface described at the MM level.
The choice of bullets are HF molecules and the target is a two dimensional BN sheet, that looks like a graphene sheet, with half of the C atoms turned into N and the other half into B atoms. In a BN sheet the atoms have of course a small charge, which we pre calculate with a QM engine (DFTB).
It is important to understand the role of bonds in this example, because the number of atoms is not constant during the simulations, as bullets are fired (and hence appear), ricochet off the surface and hence disappear after a while. They may as well stick to, or penetrate into the surface, but this is beyond the hybrid engine concept.
In this example there are fixed bond orders withing the target and within the bullets. This is because we specify GuessBonds in the two system blocks (target and bullet). When a bullet is added its bonds are automatically added. The hybrid engine itself will never guess bonds and always use what is specified on input. No bonds are ever formed between the bullet and the surface (the QM and MM regions).
#!/bin/sh
# In this example we use the hybrid engine in a molecular gun MD application, shooting HF molecules at a BN surface
# The BN slab represents the MM region and the "bullets" are the QM region
# The regions are defined in the xyz files using end of line strings (atom attributes)
# First we do two dftb calculations to get a guess of the charges to be used by the force field.
STRUCTDIR=$AMSHOME/examples/Hybrid/HybridGun/molecules
export AMS_JOBNAME=BNSlab.dftb
rm -rf $AMS_JOBNAME.results
$AMSBIN/ams << eor
Task SinglePoint
System
   GeometryFile $STRUCTDIR/BNSlab.xyz
End
Engine DFTB
EndEngine
eor
export AMS_JOBNAME=HF.dftb
rm -rf $AMS_JOBNAME.results
$AMSBIN/ams << eor
Task SinglePoint
System
   GeometryFile $STRUCTDIR/HF.xyz
End
Engine DFTB
EndEngine
eor
# now we can run our MD simulation using both mechanical and electrostatic embedding
for embedding in mechanical electrostatic
do
# because electrostatic embedding is more expensive we limit here the number of steps
steps=1400
if [ $embedding = electrostatic ]
then
  steps=300
fi
export AMS_JOBNAME=SinkBox.embedding=$embedding
rm -rf $AMS_JOBNAME.results
$AMSBIN/ams << eor
Task MolecularDynamics
System
    GeometryFile $STRUCTDIR/BNSlab.xyz
    GuessBonds true
    LoadForceFieldCharges file=BNSlab.dftb.results
End
System H2
   GeometryFile $STRUCTDIR/HF.xyz
   GuessBonds true
   LoadForceFieldCharges file=HF.dftb.results
End
RNGSeed -1341016088 83513668 1764626453 -87803069 -1149690266 1963370818 -1393571175 1985130742
MolecularDynamics
    NSteps $steps
    Trajectory
        SamplingFreq 20
    End
    InitialVelocities
        Temperature 300
    End
    AddMolecules
        System H2
        Frequency 159
        CoordsBox 0 3 0 8.57 6 7
        VelocityDirection 0.45752820 0 -0.5540656
        Velocity 0.07
        Rotate Yes
        MinDistance 3.0
    End
    Preserve
        Momentum No
        AngularMomentum No
    End
    RemoveMolecules
      Formula *
      Frequency 101
      SinkBox FractionalCoordsBox="0 1 0 1 8 1000"
   End
End
Constraints
  Atom 1
End
Engine Hybrid
   QMMM  qmRegion=qm mmEngineID=ForceField qmEngineID=dftb embedding=$embedding
   Engine dftb
   EndEngine
   Engine ForceField
      NonBondedCutoff 50 [Bohr]
   EndEngine
EndEngine
eor
done