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

Note that this example makes use of the python scripting tools PLAMS
and PISA, available in the AMS Python Stack.

"Download HybridGun.py"

   #!/usr/bin/env amspython

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

   from scm.base import ChemicalSystem, Units
   from scm.input_classes import drivers, engines
   from scm.plams import AMSJob, JobRunner, config, Settings
   import os

   # Load chemical systems from xyz files
   struct_dir = os.path.join(os.getenv("AMSHOME"), "examples", "Hybrid", "HybridGun", "molecules")
   bn_slab = ChemicalSystem.from_xyz(os.path.join(struct_dir, "BNSlab.xyz"))
   hf = ChemicalSystem.from_xyz(os.path.join(struct_dir, "HF.xyz"))
   systems = {"bn_slab": bn_slab, "hf": hf}

   # First we do two dftb calculations to get a guess of the charges to be used by the force field.
   # Setup jobs and run in parallel
   config.default_jobrunner = JobRunner(parallel=True)
   settings = Settings()
   settings.runscript.nproc = 1

   driver_dftb = drivers.AMS()
   driver_dftb.Task = "SinglePoint"
   driver_dftb.Engine = engines.DFTB()

   results_dftb = {}
   for name, system in systems.items():
       job_dftb = AMSJob(molecule=system, settings=settings + Settings({"input": driver_dftb}), name=name)
       results_dftb[name] = job_dftb.run()

   # Load forcefield charges from the dftb calculations onto the chemical systems
   for name, system in systems.items():
       results = results_dftb[name]
       charges = results.get_charges()
       system.enable_atom_attributes("forcefield")
       for atom, charge in zip(system.atoms, charges):
           atom.forcefield.charge = charge
       system.guess_bonds()
       print(f"Total Energy (hartree) {results.get_energy('hartree', 'dftb'):.8f}")
       print(f"Total Energy (eV) {results.get_energy('eV', 'dftb'):.8f}")

   # Now we can run our MD simulation using both mechanical and electrostatic embedding
   # (because electrostatic embedding is more expensive we limit here the number of steps)
   results_md = {}
   for embedding, steps in [("mechanical", 1400), ("electrostatic", 300)]:

       driver_md = drivers.AMS()
       driver_md.Task = "MolecularDynamics"
       driver_md.RNGSeed = [-1341016088, 83513668, 1764626453, -87803069, -1149690266, 1963370818, -1393571175, 1985130742]
       driver_md.Constraints.Atom = 1

       driver_md.MolecularDynamics.NSteps = steps
       driver_md.MolecularDynamics.Trajectory.SamplingFreq = 20
       driver_md.MolecularDynamics.InitialVelocities.Temperature = 300
       driver_md.MolecularDynamics.Preserve.Momentum = "No"
       driver_md.MolecularDynamics.Preserve.AngularMomentum = "No"
       driver_md.MolecularDynamics.AddMolecules.System = "HF"
       driver_md.MolecularDynamics.AddMolecules.Frequency = 159
       driver_md.MolecularDynamics.AddMolecules.CoordsBox = [0, 3, 0, 8.57, 6, 7]
       driver_md.MolecularDynamics.AddMolecules.VelocityDirection = [0.45752820, 0, -0.5540656]
       driver_md.MolecularDynamics.AddMolecules.Velocity = 0.07
       driver_md.MolecularDynamics.AddMolecules.Rotate = "Yes"
       driver_md.MolecularDynamics.AddMolecules.MinDistance = 3.0
       driver_md.MolecularDynamics.RemoveMolecules.Formula = "*"
       driver_md.MolecularDynamics.RemoveMolecules.Frequency = 101
       driver_md.MolecularDynamics.RemoveMolecules.SinkBox.FractionalCoordsBox = [0, 1, 0, 1, 8, 1000]

       driver_md.Engine = engines.Hybrid()
       driver_md.Engine.QMMM.QMRegion = "QM"
       driver_md.Engine.QMMM.QMEngineID = "DFTB"
       driver_md.Engine.QMMM.MMEngineID = "ForceField"
       driver_md.Engine.QMMM.Embedding = embedding

       driver_md.Engine.Engine[0] = engines.DFTB()
       driver_md.Engine.Engine[1] = engines.ForceField()
       driver_md.Engine.Engine[1].NonBondedCutoff = 50 * Units.conversion_factor("bohr", "angstrom")

       job_md = AMSJob(
           molecule={"": systems["bn_slab"], "HF": systems["hf"]},
           settings=settings + Settings({"input": driver_md}),
           name=f"SinkBox.embedding={embedding}",
       )
       results_md[embedding] = job_md.run()

   for result in results_md.values():
       res = result.read_rkf_section("MDResults")
       header = "{:<20} {:>16} {:>16} {:>16} {:>16}".format("Measure", "Min", "Max", "Mean", "StdDev")
       print(header)
       print("-" * len(header))
       print(
           f"{'Potential Energy':<20} {res['MinPotentialEnergy']:>16.8f} {res['MaxPotentialEnergy']:>16.8f} {res['MeanPotentialEnergy']:>16.8f} {res['StdDevPotentialEnergy']:>16.8f}"
       )
       print(
           f"{'Kinetic Energy':<20} {res['MinKineticEnergy']:>16.8f} {res['MaxKineticEnergy']:>16.8f} {res['MeanKineticEnergy']:>16.8f} {res['StdDevKineticEnergy']:>16.8f}"
       )
       print(
           f"{'Total Energy':<20} {res['MinTotalEnergy']:>16.8f} {res['MaxTotalEnergy']:>16.8f} {res['MeanTotalEnergy']:>16.8f} {res['StdDevTotalEnergy']:>16.8f}"
       )
       print(
           f"{'Temperature':<20} {res['MinTemperature']:>16.8f} {res['MaxTemperature']:>16.8f} {res['MeanTemperature']:>16.8f} {res['StdDevTemperature']:>16.8f}"
       )
