Example: Using capping atoms in a periodic system with charges
**************************************************************

Now we look at a BN chain and use charges for the MM calculation.

[image]

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.042      0.034
   B-N             1.431      1.499      0.068

   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.844721    0.004   -0.007   -0.044
   mechanical             fractional    -6.745243    0.012   -0.007   -0.042
   electrostatic               fixed    -6.749110    0.002   -0.004   -0.036
   electrostatic          fractional    -6.652361    0.010   -0.004   -0.026

   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

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

"Download PeriodicCappingWithCharge.py"

   #!/usr/bin/env amspython

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

   # Load molecule and run initial geometry optimization
   config.default_jobrunner = JobRunner(parallel=True)
   settings = Settings()
   settings.runscript.nproc = 1

   struct_dir = os.path.join(os.getenv("AMSHOME"), "examples", "Hybrid", "PeriodicCapping", "systems")

   mol = ChemicalSystem.from_xyz(os.path.join(struct_dir, "var5.xyz"))
   mol.guess_bonds()

   driver_geom_opt = drivers.AMS()
   driver_geom_opt.Task = "GeometryOptimization"
   driver_geom_opt.Engine = engines.DFTB()
   driver_geom_opt.GeometryOptimization.OptimizeLattice = True
   driver_geom_opt.GeometryOptimization.Method = "FIRE"
   driver_geom_opt.GeometryOptimization.MaxIterations = 300

   job_dftb = AMSJob(molecule=mol, settings=settings + Settings({"input": driver_geom_opt}), name="reference")
   results_dftb = job_dftb.run()

   # Get charges and add them to molecule
   mol.enable_atom_attributes("forcefield")
   charges = results_dftb.get_charges()
   for atom, charge in zip(mol.atoms, charges):
       atom.forcefield.charge = charge

   # Run cheap calculation using forcefield engine
   driver_geom_opt.Engine = engines.ForceField()
   driver_geom_opt.Engine.NonBondedCutoff = 50 * Units.conversion_factor("bohr", "angstrom")
   job_ff = AMSJob(molecule=mol, settings=settings + Settings({"input": driver_geom_opt}), name="cheap")
   results_ff = job_ff.run()
   mol = results_ff.get_main_system()

   print("begin report")
   print(
       """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:
   """
   )

   header = "{:<10} {:>10} {:>10} {:>10}".format("Distance", "QM", "MM", "Err(MM)")
   print(header)
   print("-" * len(header))

   mol_qm = results_dftb.get_main_system()
   mol_mm = results_ff.get_main_system()
   for bond, i, j in [("B-H", 0, 1), ("N-H", 2, 3), ("B-N", 0, 2)]:
       dist_qm = mol_qm.get_distance(i, j)
       dist_mm = mol_mm.get_distance(i, j)
       print(f"{bond:<10} {dist_qm:>10.3f} {dist_mm:>10.3f} {dist_mm - dist_qm:>10.3f}")

   # Run hybrid engine
   results = {}
   for embedding in ["mechanical", "electrostatic"]:
       for capping in ["fixed", "fractional"]:
           driver = copy.deepcopy(driver_geom_opt)
           driver.Properties.Gradients = True
           driver.GeometryOptimization.MaxIterations = 100

           driver.Engine = engines.Hybrid()
           driver.Engine.Capping.AllowHighBondOrders = True
           driver.Engine.Capping.Option = capping
           driver.Engine.QMMM.QMRegion = "qm"
           driver.Engine.QMMM.QMEngineID = "DFTB"
           driver.Engine.QMMM.MMEngineID = "ForceField"
           driver.Engine.QMMM.Embedding = embedding

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

           job = AMSJob(
               molecule=mol,
               settings=settings + Settings({"input": driver}),
               name=f"system.embedding={embedding}.capping={capping}.go",
           )

           results[(embedding, capping)] = job.run()

   hybrid_results = {k: (v.get_energy(engine="hybrid"), v.get_main_system()) for k, v in results.items()}

   print(
       """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:
   """
   )

   header = "{:<16} {:>16} {:>12} {:>8} {:>8} {:>8}".format("Embedding", "Capping", "Energy", "B-H", "N-H", "B-N")
   print(header)
   print("-" * len(header))

   for (embedding, capping), (energy, mol_hybrid) in hybrid_results.items():
       errs = [mol_hybrid.get_distance(i, j) - mol_qm.get_distance(i, j) for i, j in [(0, 1), (2, 3), (0, 2)]]
       print(f"{embedding:<16} {capping:>16} {energy:>12.6f} {errs[0]:>8.3f} {errs[1]:>8.3f} {errs[2]:>8.3f}")

   print(
       """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
   """
   )
   print("end report")
