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 amin=0 amax=1 bmin=0 bmax=1 cmin=8 cmax=1000
End
End
Constraints
Atom 1
End
Engine Hybrid
QMMM qmRegion=qm mmEngineID=ForceField qmEngineID=dftb embedding=$embedding
Engine dftb
EndEngine
Engine ForceField
EndEngine
EndEngine
eor
done